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I.  INTRODUCTION  AND  SUMMARY 

1 . 1  BACKGROUND 

The  objective  of  the  Systems,  Science  and  Software  ( S 3 ) 
research  program  is  to  extend  our  present  understanding  of  the 
excitation  of  seismic  waves  by  underground  explosions  and 
earthquakes.  Toward  this  objective,  we  are  conducting  theo¬ 
retical  and  empirical  studies  of  ground  motions  from  the  two 
classes  of  sources.  In  particular,  our  efforts  are  directed 
toward  the  development  of  improved  methods  for  discriminating 
between  the  seismic  signals  from  earthquakes  and  explosions 
and  the  development  of  improved  methods  for  estimating  ex¬ 
plosion  yield. 

This  report  summarizes  the  work  done  during  the  first 
three-month  period  of  the  contract. 

1.2  SUMMARY  OF  RESEARCH  DURING  THIS  QUARTER 

Our  work  during  this  quarter  has  included  research  in 
a  number  of  areas.  Research  projects  are  briefly  summarized 
below. 

Source  Calculations 
A.  Explosion  Source  Modeling 

Spherically  symmetric  one-dimensional  calculations  for 
the  SALMON  event  in  the  Tatum  salt  dome  in  Louisiana  have 
been  reported  earlier.  These  calculations,  based  on  labora¬ 
tory  material  properties  for  salt,  gave  a  reduced  displace¬ 
ment  potential  (RDP)  approximately  a  factor  of  two  lower  than 
the  RDP  based  on  measured  acceleration  and  velocity  data.  In 
the  present  reporting  period,  we  have  been  proceeding  with  a 
program  aimed  at  resolving  the  discrepancies  between  our  cal¬ 
culations  and  the  data. 


An  axi symmetric  two-dimensional  calculation  in  which 
the  material  is  initially  assumed  to  be  in  uniaxial  strain 
(including  a  shear  stress)  has  been  completed.  This  prestress 
condition  gives  an  asymmetric  radiation  pattern,  with  a  30 
percent  increase  in  static  displacement  at  measuring  stations 
at  shot  level,  but  a  60  percent  decrease  vertically.  Further 
investigation  has  linked  this  radiation  pattern  directly  to 
the  in  situ  shear  stress.  Another  two-dimensional  calculation 
has  been  completed  which  included  the  anhydrite,  limestone, 
and  sediments  above  the  working  point.  We  concluded  that 
these  layers  could  not  significantly  increase  the  horizontal 
static  displacement  at  shot  level. 

An  intensive  literature  review  of  the  material  proper¬ 
ties  of  salt  rocks  is  near  completion.  We  have  found  ex¬ 
tensive  in  situ  and  laboratory  data  for  salt  rocks,  much  of 
it  developed  in  connection  with  the  Waste  Isolation  Pilot 
Plant  (WIPP) .  Based  on  this  data,  we  have  introduced  thermal 
softening  into  our  salt  constitutive  model.  Calculations  have 
been  made  using  thermal  softening  plus  a  more  judicious  choice 
of  elastic  constants  based  on  in  situ  seismic  data,  and  a 
failure  surface  which  envelopes  the  great  mass  of  laboratory 
salt  data.  The  latest  one-dimensional  calculation  gives  an 
RDP  only  20  percent  lower  than  the  measured  data  and  peak 
stresses  and  velocities,  well  within  the  scatter  of  the  mea¬ 
surements.  Thus,  it  appears  that  the  strain  rate  dependent 
constitutive  behavior  of  salt  should  be  included  in  our  model. 

B.  Earthquake  Modeling  on  the  ILLIAC  IV  Computer 

This  work  is  summarized  in  Section  1.6  and  described 
in  detail  in  Section  V. 


Discrimination 

This  work  is  summarized  in  Section  1.7  and  described 


in  detail  in  Section  VI. 


Yield  Determination 


V# 


A.  Time  Domain  m^ 

Some  procedures  for  improving  the  accuracy  and  conve¬ 
nience  of  the  measurement  of  time  domain  m^  are  proposed  in 
summary  form  in  Section  1.3.  A  detailed  discussion  is  given 
in  Section  II. 

“**  *•*•••* — wMCi _  «  —  - '  •  — . . .  — - - - 

A 

B.  mjj  An  Automated  Spectral  Magnitude 

A  new  body  wave  magnitude,  m^,  is  proposed  in  S'-  ion 
III.  This  section  is  summarized  in  Section  1.4. 

C.  mb  for  Synthetic  Seismograms 

This  work  is  summarized  in  Section  1.5  and  is  dt  .bed 
in  detail  in  Section  IV. 

Small-Scale  Experiments 

We  have  now  successfully  test  fired  two  of  the  minia¬ 
ture  spherical  charges  for  use  in  the  model  experiments.  The 
first  spheres  constructed  would  not  detonate  because  of 
impurities  in  the  explosive,  PETN;  this  problem  has  now  been 
corrected.  A  repeatable  set  of  construction  procedures  is 
being  developed.  The  detailed  design  of  the  initial  proof 
experiment  has  been  started;  the  objective  of  this  experiment 
will  be  duplication  of  the  results  obtained  in  the  1977  tests. 

Selected  Geological  Studies 

A  meeting  of  the  Russian  Geophysics  Study  Group  was 
held  at  SDAC  on  December  11  and  12.  The  group  included 
representatives  of:  (1)  government  agencies  who  are  funding 
seismic  discrimination  research  (ARPA,  AFTAC,  CIA) ,  (2)  con¬ 
tractors  who  are  carrying  out  research  in  seismic  discrimina¬ 
tion  (S3,  Teledyne-Geotech),  (3)  groups  that  are  monitoring 
the  Russian  earth  science  literature  (USGS,  RAND,  LLL) ,  and 


3 


b 
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(4)  geophysical  consultants  who  are  knowledgeable  about  the 
Russian  earth  science  research  in  their  areas  of  expertise 
(Alan  Ryall,  Shelton  Alexander,  Rob  Wesson,  George  Keller, 
George  Woolard,  Alex  Malahoff) .  The  general  objective  of  the 
meeting  was  to  assess  the  adequacy  of  the  existing  geophysical 
data  base  to  be  used  in  seismic  calculations  for  potential 
Russian  testing  areas.  The  opening  session  centered  on  a 
series  of  presentations  by  the  contractors  in  which  they 
attempted  to  provide  quantitative  estimates  of  the  sensitivity 
of  the  calculations  to  uncertainties  in  the  specification  of 
the  geophysical  parameters  of  the  source  region  and  propaga¬ 
tion  paths.  This  was  followed  by  a  series  of  presentations 
by  groups  that  are  monitoring  the  Russian  earth  science 
research  in  which  an  attempt  was  made  to  summarize  the  status 
of  various  efforts  to  define  values  for  the  physical  param¬ 
eters  of  interest.  This  was  followed  by  an  open  discussion 
in  which  gaps  in  the  existing  data  base  were  identified  and 
possible  means  of  filling  these  gaps  were  suggested  by  the 
various  consultants.  Minutes  of  the  meeting  are  currently  in 
preparation  and  should  be  available  for  distribution  in  the 
coming  month. 


Magnitude-Yield  Improvements 

Initial  contact  has  been  made  with  AFTAC  and  sample 
ground  motion  spectra  have  been  obtained.  A  quantitative 
description  of  the  data  and  a  list  of  the  events  and  stations 
for  which  data  are  available  is  now  being  compiled. 


Ground  Motion  Analysis 

A  systematic  search  for  all  published  free-field  seis¬ 
mic  data  for  events  in  tuff,  rhyolite  and  alluvium  is  under¬ 
way  . 
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1.3  SUMMARY  OF  SECTION  II:  "SEISMIC  WAVEFORMS  AND  TIME 

DOMAIN  mb" 

Body  wave  magnitude,  mb,  is  an  important  single  param¬ 
eter  used  to  describe  a  large  body  of  data,  seismograms  from 
many  stations.  Conventional  is  based  on  direct  measure¬ 
ments  made  by  an  analyst  from  analog  playouts  of  the  data  and 
includes  a  correction  for  the  response  of  the  seismometer  at 
the  apparent  period  of  the  phase  measured.  When  digital  data 
are  available,  as  is  increasingly  the  case,  this  procedure 
is  unnecessarily  cumbersome  and  prone  to  error. 

We  suggest  a  semi-automated  procedure  that  essentially 
eliminates  measurement  errors.  We  first  filter  all  seismograms 
so  they  appear  as  if  recorded  by  the  same  seismometer.  This 
removes  a  source  of  systematic  differences  that  can  be  0.2  m^ 
units  or  more.  The  time  and  amplitude  of  peaks  on  the  record 
are  then  determined  automatically  by  a  parabolic  fit  to  a 
moving  three  point  window.  The  phase  to  be  used  for  mb  is 
selected  by  direct  examination  of  the  waveform  and  its  period 
and  amplitude  are  read  from  a  table. 

The  semi -automated  procedure  outlined  above  is  demon¬ 
strated  by  applications  to  HNME  recordings  of  eleven  Pahute 
Mesa  explosions.  The  resulting  m^  are  compared  to  those  given 
by  SDAC .  A  more  complete  description  of  this  work  is  given 
by  Bache  (1979) . 

1.4  SUMMARY  OF  SECTION  III:  "mb ,  AN  AUTOMATED  SPECTRAL 

MAGNITUDE" 

In  this  section  our  objective  is  to  determine  and  test 

A 

an  automated  spectral  magnitude  which  we  call  m^.  This  mag¬ 
nitude  is  based  on  the  spectral  amplitude  in  a  narrow  window 
in  both  time  and  frequency.  The  concept  of  the  m^  has  emerged 
from  parallel  development  of  an  automated  algorithm  for  nuclear 
explosion  detection  and  discrimination.  The  basic  procedure 


on  which  we  base  the  detection,  discrimination  and  mb  deter¬ 
mination  is  the  same  and  is  contained  in  the  MARS  (Multiple 
Arrival  Recognition  System)  program. 

Like  any  magnitude  measure,  the  mb  is  an  empirical 
measure  that  must  be  developed  by  testing  on  actual  data.  The 
data  set  used  in  Section  III  includes  HNME  recordings  of 
eleven  events.  We  describe  the  application  of  MARS  to  these 

A 

data  in  considerable  detail  and  show  how  the  mb  is  computed. 

The  basic  operation  is  filtering  with  a  Gaussion  filter  that 
is  narrow  in  both  time  and  frequency.  Empirical  tests  lead 
to  a  preferred  filter  which  has  a  width  of  0.16  Hz  or  4 
seconds.  Here  the  width  is  defined  to  include  the  central 
region  that  contains  68  percent  of  the  area  under  the  Gaussion 
filter. 

Bache  (1979)  discussed  the  mb  for  a  much  larger  data 
sample  including  seismograms  from  six  stations.  He  concludes 
that  mb  is  at  least  as  good  a  magnitude  measure  as  the  most 
carefully  determined  time  domain  mb,  at  least  for  high 
signal/noise  data. 

1.5  SUMMARY  OF  SECTION  IV:  "mb  FOR  SYNTHETIC  SEISMOGRAMS" 

To  improve  understanding  of  the  mb,  we  apply  it  to 
synthetic  seismograms  that  closely  match  the  HNME  recordings 
of  three  Pahute  Mesa  explosions.  The  construction  of  these 
synthetic  seismograms  has  considerable  intrinsic  interest 
and  this  is  described. 

1.6  SUMMARY  OF  SECTION  V:  "THREE-DIMENSIONAL  EARTHQUAKE 

MODELING " 

In  Section  V,  we  summarize  our  earthquake  source  modeling 
research  using  the  ILLIAC  computer.  Much  of  this  work  was 
accomplished  under  a  separate  contract  with  Air  Force  Geophysics 
Laboratory  and  is  reported  in  detail  by  Day,  et  al . ,  (1978). 

The  results  presented  in  this  section  are  a  summary  of  those 


given  in  our  report  to  AFGL.  We  also  discuss  our  plans  for 
research  with  this  model  under  this  contract. 

Two  three-dimensional  finite  difference  calculations 
were  performed  for  faults  with  constant  rupture  velocity  in 
a  uniform  prestress  field.  In  one  case,  the  medium  was 
modeled  as  perfectly  elastic;  in  the  second  case,  the  medium 
was  elastic-plastic.  Accuracy  of  the  numerical  method  has 
been  verified  by  comparison  to  an  exact,  closed-form  solution. 
Near-  and  far-field  seismic  signals  have  been  obtained  from 
the  numerical  solutions.  Comparing  the  radiated  field  for  the 
finite  difference  solution  to  that  for  the  source  model  of 
Archambeau  has  helped  clarify  the  physical  interpretation  of 
the  parameters  of  the  Archambeau  spherical  source  model. 

The  inclusion  of  plastic  yielding  in  the  second  finite 
difference  calculation  is  a  significant  step  toward  incorpo¬ 
rating  realistic  rock  mechanics  into  the  fault  model.  It  was 
found,  however,  that  this  simple  form  of  plasticity  did  not 
reduce  the  large  velocity  peaks  observed  on  the  fault  plane 
in  the  elastic  case.  This  result  supports  previous  studies 
which  have  suggested  that  an  abrupt  stress  drop  is  incon¬ 
sistent  with  bounded  velocities. 

1.7  SUMMARY  OF  SECTION  VI:  "DISCRIMINATION  EXPERIMENT" 

Two  modifications  were  made  to  the  MARS  program  that 
we  are  using  in  the  discrimination  experiment.  One  of  these 
involved  the  development  of  an  algorithm  for  estimating  and 
correcting  for  the  effects  of  local  seismic  noise  (i.e., 
either  background  noise  or  signal  coda)  on  a  transient  sig¬ 
nal.  The  second  modification  consists  of  averaging  weighted 
magnitude  estimates  at  several  different  frequencies  over  a 
high  and  low  frequency  band.  The  weights  are  based  on  the 
ratio  of  signal  power  to  noise  power. 

Event  seismograms  were  processed  with  the  modified 
MARS  program  for  classified  stations.  Preliminary  results 


based  on  52  events  recorded  at  one  or  more  of  these  stations 
indicate  that  the  variable  frequency  magnitude  (VFM)  approach 
can  successfully  discriminate  events  down  to  small  magnitude 
levels  at  stations  with  low  background  noise  levels.  There 
is  also  an  indication  that  stations  located  over  high-Q  upper 
mantle  regions  (e.g.,  shields)  provide  better  separation  of 
earthquakes  and  explosions  that  stations  located  in  tectonic 
regions  of  high  heat  flow  and  large  positive  (slow)  travel¬ 
time  residuals  (by  inference,  low-Q  upper  mantle) .  This 
tentative  result  will  be  examined  in  more  detail  when  the 
total  data  base  (classified  and  unclassified  stations)  is 
completed. 


II.  SEISMIC  WAVEFORMS  AND  TIME  DOMAIN 


"b 


2.1  INTRODUCTION 


In  this  section  we  present  a  method  for  determining 
time  domain  m^  which  uses  conventional  formulas,  but  takes 
advantage  of  the  availability  of  digital  data  to  reduce  mea¬ 
surement  errors.  The  method  will  be  applied  to  recordings  of 
Pahute  Mesa  explosions.  The  resulting  m^  are,  we  believe,  the 
best  time  domain  values  possible  for  each  record.  We  will  com¬ 
pare  these  to  those  given  by  SDAC  for  the  same  records.  In 
later  sections  we  will  be  discussing  a  spectral  magnitude 


called  m^.  These  measurement  error  free  m^  are  needed  for 
comparison  to  the  m^. 


The  data  to  be  analyzed  are  short  period  recordings 
of  eleven  Pahute  Mesa  events  from  the  SDAC  station  HNME.  The 
event  locations  are  shown  in  Figure  1  on  a  plot  which  includes 
the  outline  of  the  Silent  Canyon  Caldera  (Orkild,  et  al. , 
1969).  Only  STILTON  lies  outside  the  Caldera  boundary. 


Source  information  for  the  eleven  events  is  summarized 
in  Tables  1-3.  The  date,  location  and  depth  are  given  in 
Table  1.  In  Table  2  we  summarize  the  known  elastic  proper¬ 
ties  of  the  source  material  and  the  overburden.  In  Table  3 
are  three  estimates  of  the  source-to-surf ace  travel  time 
which  is  roughly  half  the  delay  time  between  the  teleseismic 
P  and  pP.  One  estimate  is  the  quotient  of  the  depth  and  the 
sonic  velocity  of  the  overburden.  For  the  second  estimate 
we  use  the  mean  overburden  velocity  (oi’OVB)  for  the  eleven 
events  rather  than  the  specific  value  for  each  event.  The 
third  estimate,  the  "observed",  is  based  on  the  measured  ar¬ 
rival  time  on  the  accelerogram  at  surface  ground  zero.  The 
references  for  the  accelerogram  records  are  the  L  ,  LASL  or 
Sandia  memoranda  indicated  in  the  tables. 


We  will  be  concerned  with  the  conventional  magnitude, 
m^,  for  these  events.  For  each  recording  this  is  determined 
according  to 
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=  log  |  +  B  (A ) 

where  A  and  T  are  the  amplitude  and  period  for  a  particular 

phase  on  the  seismogram  with  A  corrected  for  the  instrument 

amplification  at  the  period  T.  The  B(A)  is  the  standard 

Gutenberg-Richter  distance  correction  which  is  3.23  for  the 

NTS-HNME  range  of  36.7  degrees.  We  will  make  no  correction 

for  distance  variations  between  events  since  these  are  less 

than  +  0.01  m^  units  in  most  cases.  The  Sc  represents  source- 

path-station  corrections  which  may  be  determined  empirically 

and  included.  However,  we  will  take  S  to  be  zero  in  this 

c 

case. 

The  instrument  amplification  can  be  an  important  fac¬ 
tor  in  the  m^  formula.  These  eleven  events  were  recorded  by 
two  different  seismometers.  The  1975  events  were  recorded  by 
the  18300  while  the  KS36000  was  operating  for  the  others. 

The  amplification  curves  for  these  two  seismometers  are 
plotted  in  Figure  2.  The  KS36000  amplitude  response  was  pro¬ 
vided  tc  S3  in  a  13  March  1978  memorandum  by  Captain  M.  J. 
Shore,  VSC.  The  phase  response  was  obtained  by  assuming  the 
instrument  was  a  minimum  phase  system.  The  amplitude  and 
phase  are  then  related  through  Hilbert  transforms.  The  scheme 
for  computing  the  phase  response  from  the  amplitude  response 
is  that  of  Bolduc,  et  al.  (1972).  For  the  18300  ,  amplitude 
and  phase  response  data  for  HNME  were  provided  in  a  letter 
by  R.  W.  Alewine,  VSC.  The  Hilbert  transform  program  was 
used  to  recompute  the  phase  which  turned  out  to  be  a  smoothed 
vision  of  the  given  data.  The  HNME  station  logs  were  checked 
against  the  nominal  18300  and  KS36000  curves  and  were  found 
to  be  in  good  agreement  for  all  eleven  events. 

In  the  remainder  of  this  section  we  will  present  the 
waveforms  as  they  were  given  to  us  and  list  the  m^  given  by 
the  SDAC  event  reports.  These  values  will  be  compared  to 


RELATIVE  AMPLIFICATION 


PERIOD  (sec) 


Figure  2.  Relative  amplitude  response  for  two  instruments 
used  at  the  HNME  site. 
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time  domain  m^  computed  by  us  using  Eq.  (1)  ,  but  with  a  new 
method  for  determining  A  and  T.  We  call  this  a  semi-auto¬ 
mated  time  domain  m^. 

2.2  A  SEMI -AUTOMATED  TIME  DOMAIN 

To  compute  m^  from  Eq .  (1)  we  need  to  measure  A  and  T 
and  this  is  usually  done  by  hand  from  analog  playouts  of  the 
data.  The  most  important  errors  inherent  in  this  procedure 
may  be  divided  into  two  categories : 

•  It  is  difficult  to  measure  the  period  with  an 
error  much  less  than  0.1  to  0.2  seconds.  The 
instrument  response  is  rapidly  changing  in  the 
0.7  to  1.4  second  range  where  the  measurements 
are  usually  made  (Figure  2) .  Fortunately,  the 
approximately  inverse  proportionality  between  T 
and  the  amplification  does  tend  to  compensate  for 
errors  in  T  and  reduce  their  effect  in  many  cases. 

•  Even  if  we  could  measure  the  period  precisely,  the 
signal  is  not  monochromatic.  Thus,  correcting  to 
true  ground  motion  by  dividing  by  the  instrument 
amplification  at  the  period  T  must  introduce  error. 
For  a  given  waveform  the  error  is  constant,  so  for 
similar  signals  recorded  by  identical  instruments, 
the  main  effect  is  a  baseline  shift.  However, 

the  error  is  compounded  when  data  from  instruments 
with  different  response  curves  are  mixed. 

We  wanted  to  compute  time  domain  m^  values  with  the 
errors  minimized  to  the  extent  possible.  This  has,  of  course 
a  considerable  amount  of  inherent  interest.  Our  original 
motivation  was,  however,  to  remove  as  many  obscuring  effects 
as  we  could  for  comparison  of  time  domain  m^  with  our  new  m^. 
Also,  as  long  as  digital  data  were  available,  we  could  see 
little  point  in  making  measurements  from  analog  playouts. 


Our  semi- automated  time  domain  m^  involves  two  steps 
which  may  or  may  not  be  used  together.  First,  we  use  the 
computer  to  make  the  measurements.  The  data  are  passed 
through  a  sliding  three-point  window.  Whenever  the  center 
point  is  greater  than  or  less  than  both  end  points ,  the  three 
points  are  fit  by  a  parabola.  The  time  and  amplitude  of  the 
peak  are  saved.  The  period,  T,  is  computed  from  the  time 
difference  between  adjacent  peaks.  Using  the  apparent  instru¬ 
ment  response,  log  (A/T)  is  then  computed  for  each  cycle  on 
the  record.  To  obtain  an  m^ ,  the  analyst  need  only  inspect 
the  waveform  to  select  a  phase  and  then  read  log  A/T  from  a 
printed  table.  This  algorithm  is  contained  in  a  subroutine 
called  PFIT  which  is  routinely  applied  to  synthetic  or  ob¬ 
served  data  for  which  m^  or  Ms  information  is  desired. 

We  have  outlined  a  method  to  make  time  domain  measure¬ 
ments  convenient  and  accurate  when  digital  data  are  available. 
What  about  the  correction  for  instrument  response?  We  be¬ 
lieve  the  effect  of  this  correction  is  most  consistent  when 
all  recordings  are  prefiltered  to  appear  as  if  recorded  by 
the  same  instrument.  That  is,  we  Fourier  transform  the  sig¬ 
nal,  multiply  by  the  quotient  of  the  actual  and  standard  re¬ 
sponses  and  inverse  transform.  The  new  recordings  can  then 
be  processed  by  PFIT.  In  later  sections  we  will  show  that 
changing  the  seismograms  to  a  standard  instrument  can  remove 
a  significant  source  of  error  which  is  far  from  normally 
distributed. 

2 . 3  HNME  DATA 

The  HNME  seismograms  are  shown  in  Figure  3  as  they 
were  recorded.  They  are  arranged  on  the  page  in  order  of  in¬ 
creasing  depth  as  given  in  Table  1.  As  noted  in  Section  2.1, 
these  seismograms  were  recorded  with  two  different  instru¬ 
ments.  The  main  difference  between  the  two  (Figure  2)  is  at 
short  periods  where  the  18300  is  significantly  larger. 


COLBY 


Figure  3.  The  HNME  recordings  of  eleven  Pahute  Mesa  events 
are  arranged  according  to  increasing  depth.  The 
asterisk  denotes  events  recorded  by  the  18300 
seismometer. 
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In  Figure  4  we  replot  the  seismograms  arranged  accord¬ 
ing  to  increasing  source-surface  travel  time,  as  given  in 
Table  3.  For  all  but  STILTON  and  FONTINA,  the  travel-time 
estimate  is  based  on  the  arrival  time  at  surface  ground  zero. 
The  mean  velocity  of  the  overburden  (aQYg)  was  used  to  esti¬ 
mate  the  time  for  the  other  two  events . 


In  Figure  4  the  seismograms  recorded  by  the  18300  in¬ 
strument  were  filtered  so  they  appear  as  if  they  were  recorded 
by  a  KS36000  seismometer.  The  effect  is  to  remove  some  of 
the  high  frequency  details  from  the  18300  recorded  events 
(compare  to  Figure  3)  . 

Displayed  as  in  Figure  4,  record  characteristics  re¬ 
lated  to  the  source-surface  travel  time  become  apparent. 

The  nature  and  causes  of  these  differences  are  subjects  for 
a  separate  study.  It  is  clear  that  the  first  few  cycles  of 
the  seismograms  contain  at  least  three  distinct  arrivals. 
First,  there  is  the  direct  P  wave.  This  is  followed  a  second 
or  so  later  by  an  arrival  we  would  identify  as  pP.  A  break 
in  the  waveform  associated  with  this  arrival  can  be  seen  on 
the  MUENSTER,  COLBY,  ESTUARY  and  STILTON  records,  but  it  is 
less  clear  on  the  others.  The  third  arrival  is  very  distinct 
and  has  a  clearly  apparent  depth  moveout.  It  arrives  about 
two  seconds  after  P,  too  late  to  be  pP.  Other  authors  have 
noted  the  presence  of  a  secondary  phase  like  this  on  similar 
explosion  records  and  have  associated  it  with  spall  closure 
(e.g„,  Springer  (1974)).  Following  Springer,  we  will  denote 
this  phase  by  P  . 

In  Table  4  we  list  two  estimates  for  for  the 
KS36000  recorded  events  and  three  estimates  for  the  other 
five  events.  The  SDAC  data  are  taken  directly  from  the  SDAC 
event  reports. 

The  S3  measurements  were  made  using  the  procedure 
described  in  the  preceding  section.  The  PFIT  routine  was 
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10  Seconds 

The  HNME  recordings  are  arranged  according  to  in¬ 
creasing  source-surface  travel  time  (At  from  Table 
3)  .  All  18300.  recorded  events  (marked  with  an 
asterisk)  have  been  filtered  to  appear  as  if 
recorded  by  the  KS36000  instrument. 


Figure  4. 


TABLE  4 


CONVENTIONAL  for  HNME  RECORDINGS* 

S 3  Me  as  uremen ts 

SDAC**  Data  18300  KS36000 


Event 

"b 

Period 

mb 

Period 

% 

Peri od 

STILTON 

5.55 

0.7 

5.60 

0.8 

5.86 

1.1 

POOL 

6.37 

1.3 

6.29 

1.3 

ESTUARY 

6.25 

1.  5 

6.07 

1.3 

TYBO 

5.37 

1.4 

6.20 

1.2 

6.26 

1.2 

MAST 

6.21 

1.1 

6.  ]C 

1.0 

6.25 

1.1 

CHESHIRE 

6.03 

1.0 

6.02 

1.1 

CAMEMBERT 

6.25 

1.0 

6.24 

1.0 

6.37 

1. 1 

MUENSTER 

6.39 

0.8 

6.49 

0.9 

COLBY 

6.38 

0.9 

6.50 

1.3 

KASSERI 

6.46 

1.0 

6.50 

1.1 

6.59 

1.2 

FONTINA 

6.48 

1.3 

6.43 

1.  3 

★ 

All  period  measurements  were  made  from  the  first  peak  to  the 
second  peak  as  shown  above  the  TYBO  record  in  Figure  4.  The 
amplitudes  were  measured  from  first  trough  to  first  peak. 

*  * 

From  SDAC  event  reports  —  authors:  J.  R.  Woolson,  K.  J. 
Hill,  D.  D.  Solari,  M.  S.  Dawkins,  M.  D.  Gillespie,  R.  R. 
Baums tark,  R.  J.  Markle,  D»  J.  Reinbold. 


21 


applied  to  the  records  as  they  were  recorded  (Figure  3)  and 
after  filtering  so  all  include  the  KS36000  response.  For  the 
five  events  with  m^  from  both  instruments,  the  differences 
are  striking.  The  T  from  the  18300  recordings  are  0.14 
seconds  shorter,  on  the  average,  than  the  T  from  the  KS36000 
records.  As  a  result,  the  m^  are  an  average  of  0.14  units 
smaller.  The  differences  are  greatest  for  STILTON  and  the 
reason  can  easily  be  seen  by  comparing  the  two  waveforms 
(Figures  3  and  4) .  However,  the  averages  for  the  other  four 
events  are  0.10  seconds  and  0.11  m^  units,  still  quite 
large. 

The  differences  between  our  (mixed  instrument)  and 
those  given  by  the  SDAC  reports  can  mostly  be  explained  by 
differences  in  the  period.  Our  period  measurements  are  very 
accurate  since  they  were  done  automatically  by  PFIT.  A  major 
difference  occurs  for  TYBO  where  the  amplitude  in  the  SDAC 
report  must  be  in  error.  Using  copies  of  the  station  logs  and 
digital  playouts  of  the  calibration  steps,  we  recalibrated 
all  the  data.  Thus,  the  gain  we  used  is  probably  not  identi¬ 
cal  to  that  used  by  SDAC.  Ignoring  TYBO,  the  differences 
between  the  SDAC  and  ours  obtained  from  recalibrated  data 
with  PFJT  range  from  -0.18  to  0.12  m,  units. 


cm 


III.  m^,  AN  AUTOMATED  SPECTRAL  MAGNITUDE 
3.1  INTRODUCTION 

In  this  section  we  describe  a  new  magnitude  measure 

A 

that  we  call  m^.  This  magnitude  is  based  on  the  spectral  amp¬ 
litude  in  a  narrow  window  in  both  time  and  frequency.  Our 
intention  is  for  the  m^  to  provide  an  estimate  of  the  spectral 
energy  of  the  direct  P  wave  arrival.  This  is  precisely  what 
is  needed  to  infer  explosion  yield.  Empirical  tests  suggest 
that  the  m^  is  a  stable  measure  of  signal  energy  and  should  re¬ 
duce  the  scatter  in  network  m^  estimates.  Another  important 
advantage  is  that  the  m^  can  be  computed  as  part  of  an  auto¬ 
mated  seismic  data  analysis  system. 

An  important  problem  for  VSC  is  the  development  of  re¬ 
lationships  between  body  wave  magnitude  and  nuclear  explosion 
yield.  Empirical  work  directed  toward  this  objective  has 
focussed  on  time  domain  m^  measured  from  the  analog  recordings 
(e.g.,  Dahlman  and  Israelson,  1977).  The  idea  of  replacing  such 
time  domain  measures  with  some  spectral  magnitude  has  been 
around  a  long  time.  For  both  earthquakes  and  explosions  the 
spectral  energy  in  the  direct  P  wave  is  probably  the  best  in¬ 
dicator  of  the  source  energy. 

There  are  several  reasons  why  a  spectral  magnitude 
has  not  emerged  to  replace  the  conventional  m^.  Foremost  is 
probably  the  fact  that  the  data  must  be  in  digital  form  and 
be  conveniently  accessible  for  analysis.  There  is  now  a  large 
enough  digital  data  base  that  this  should  no  longer  be  a  prob¬ 
lem.  Another  reason  is  that  a  really  satisfactory  algorithm 
for  computing  spectral  has  not  yet  emerged.  However,  prom¬ 
ising  results  have  been  obtained  in  some  recent  work  at  the 
Seismological  Institute,  Uppsala  (Shapira  and  Kulhanek,  1978) 
and  at  VSC  (V7oodward,  personal  communication)  .  In  the  former 
study  the  log  A/T  in  (1)  is  replaced  by  the  mean  of  the  Fourier 
spectrum  of  the  signal  over  a  narrow  frequency  range;  the 
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authors  suggest  0.47  to  0.78  Hz.  They  find  that  this  magni¬ 
tude  has  less  scatter  than  conventional  m^  over  the  five 
station  Swedish  seismograph  station  network.  The  spectral 
magnitude  used  by  Woodward  at  VSC  is  similar  to  that  of 
Shapira  and  Kulhanek  and  he  finds  less  scatter  in  magnitude 
versus  log  yield  for  NTS  explosions. 

The  spectral  magnitudes  mentioned  in  the  previous  para¬ 
graph  are  based  on  the  Fourier  spectrum  of  some  signal  window 
that  is  much  broader  than  the  P  wave  itself.  Later  arriving 
phases  then  can  have  a  strong  effect  on  the  spectral  amplitude. 
Knowledge  of  this  disadvantage  has  probably  dampened  enthusiasm 
for  magnitudes  based  on  spectral  information. 

The  idea  for  the  mfa  emerged  from  our  development  of 
the  automated  algorithm  for  nuclear  explosion  detection  and 
discrimination  on  which  the  MARS  (Multiple  Arrival  Recognition 
System)  program  is  based.  The  MARS  algorithm  is  based  on 
analysis  of  seismic  waveforms  with  a  suite  of  narrow-band 
filters.  For  detection  and  discrimination  the  amplitude  of 
the  filter  output  is  used  to  define  a  VFM  (variable  frequency 

/v 

magnitude) .  Our  m^  is  based  on  the  VFM. 

In  the  sequel  we  first  describe  the  MARS  algorithm 

A 

and  its  application  to  define  m^.  As  examples,  the  HNME 
seismograms  from  Section  II  are  processed.  A  clearer  under- 

A 

standing  of  the  sensitivity  of  mfa  to  the  waveform  character 
can  be  obtained  by  processing  synthetic  seismograms.  Some 
results  of  such  an  exercise  are  presented  in  Section  IV. 

3.2  MARS  ANALYSIS  AND  THE  DEFINITION  OF  VFM 

The  MARS  program  has  been  developed  at  S3  during  the 
past  several  years.  The  algorithm'  and  much  of  the  program 
design  is  due  to  Professor  C.  B.  Archambeau,  an  S3  consultant. 
The  program  development  has  been  primarily  by  Dr.  J.  F.  Masso. 


Applications  have  included  nuclear  explosion  detection  and 
discrimination  (Savino  and  Archambeau,  1974;  Savino,  et  al. , 
1974>  Savino,  et  al. ,  1975;  Rodi,  et  al . ,  1978),  decomposi¬ 
tion  of  multiple  events  (Lambert,  et  al. ,  1977;  Lambert  and 
Bache,  1977)  and  analysis  of  surface  wave  dispersion  (Bache, 
et  al. ,  1978).  A  version  of  the  program  is  currently  opera¬ 
tional  on  the  Network  Event  Processor  at  VSC.  The  MARS  algo¬ 
rithm  also  promises  to  be  a  powerful  event  detector  and,  in 
short,  a  tool  that  can  form  the  core  of  a  complete  seismic 
data  analysis  system. 

Briefly  described,  the  MARS  algorithm  is  as  follows. 

The  seismogram  is  Fourier  transformed  and  filtered  by  a  narrow- 
band,  Gaussian  filter.  The  Hilbert  transform  of  the  filtered 
spectrum  is  then  constructed  and  inverse  transformed  to  obtain 
the  envelope  function.  The  peaks  of  the  envelope  function 
indicate  phase  arrivals  of  energy  in  a  narrow-band  of  fre¬ 
quencies  near  the  filter  center  frequency.  The  arrival  times 
are  accurately  preserved  in  the  times  of  these  peaks  and  the 
relative  amplitudes  of  the  peaks  reflect  the  relative  ampli¬ 
tudes  of  arriving  pulses  of  energy  in  the  frequency  band. 

The  flow  of  operations  in  the  MARS  program  is  summarized  in 
Figure  5. 

The  output  of  the  narrow-band  filter  processing  that 

forms  the  heart  of  the  MARS  program  is  a  table  including  the 

envelope  peak  amplitudes  (A^)  and  associated  arrival  times 

(tg)  for  each  center  frequency  (f^) .  Generally  there  are 

many  A,  ,  t  pairs  for  each  filter,  depending  on  the  complexity 
k  g 

of  the  signal. 

Let  us  look  more  closely  at  the  narrow-band  filter. 

As  given  in  Figure  5,  the  functional  form  is 

F(k>  (w)  =  J2e-a(f-fk)2  (5) 
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NARROW  BAND  FILTER 


^|P 


CALCULATE  QUADRATURE  SPECTRUM 
BY  HILBERT  TRANSFORM 

Y(k)(ui)  =  -i  sgn  (w)  X(k)  (o>) 

IFFT  X  (cu)  ,  Yfu)  TO  OBTAIN 

Z(k)  (t)  *  X(k)  (t)  +  i  Y(k)  (t) 

=  A(k)  (t)  exp  (i*  (k)  (t)  ) 

FORM  ENVELOPE  FUNCTION  IN  TIME  DOMAIN 
A  {k)  (t)  =  |Z(k)(t}|  =  J  [X(k)  ( t)  ]  2  +  [Y(k)  (t)  ]2 


CALCULATE  INSTANTANEOUS  PHASE 
$  (k)  (t)  =  u,t  +  4>  (k)  (t) 


=  arctan  [Y  (k)  (t) /X (k)  ( t)  ] 


AND  FREQUENCY 


-H 


(k) 


.  d<t 

X,  -r-r 

.<  at 


(k) 


Fiaure  5. 


(continued) 


where 


a 


SLn  2  Q2 


Viewed  as  a  Gaussian  probability  distribution,  this  filter  has 
standard  deviation 


VXn2  Q 


which  means  that  68.3  percent  of  the  area  under  the  filter  is 
contained  in  the  frequency  band  -  of  £  +  of . 

It  is  easily  shown  that  the  time  domain  expression  of 
the  filter  (5)  is 


F(k)  (t)  =  cos  <27Tfkt)  e  a  (7) 

The  Gaussian  envelope  function  then  has  standard  deviation 


•JTn2  Q 


(8) 


We  see  that  the  width  of  the  narrow-band  filter  is 
controlled  by  the  parameter  Q  which  appears  in  the  definition 
of  a.  From  (6)  and  (8)  it  is  clear  that  the  narrower  the  filter 
in  frequency,  the  wider  in  time  and  vice  versa.  For  high  Q  we 
get  an  excellent  estimate  of  the  spectral  amplitude  for  a  broad 
time  window;  there  will  be  few  A,  ,  t  pairs  for  each  filter. 


For  low  Q  the  spectral  estimate  is  less  precise  but  the  t  are 
more  accurate  and  there  are  generally  many  peaks  for  each  filter. 

A  choice  of  a  preferred  Q  was  made  by  trial  and  error 
testing  of  different  values  on  both  synthetic  and  observed 
seismograms.  For  short  period  data  sampled  20  times/second, 
we  have  adopted  the  following  Q  values: 


Q  -  15  V 

f  k~0  •  35 

/ 

Q  =  5.25, 

fk<0.35 

• 

^  (k) 

The  narrow  band  filters  F  (to)  ,  for  the  Q  values 
given  in  (9)  and  the  f^  used  in  our  analysis  are  plotted  in 
Figure  6.  For  f^s0.35  we  compute  from  (6)  that 

=  0.08  Hz 

Thus,  the  width  of  the  filter  within  one  standard  deviation  is 
0.16  Hz.  For  f^<0.35  equation  (6)  gives  af  =  0.23  f^  and  the 
filter  is  narrower. 

The  corresponding  time  domain  filters  are  plotted  in 
Figure  7.  From  (8)  we  have  for  f^O.35  that 

aT  =  1.99  seconds. 

The  filter  width  is  therefore  about  four  seconds  for  these  f^, 
using  standard  deviation  as  the  criterion,  and  somewhat  wider 
at  lower  center  frequencies.  Clearly,  this  is  not  sufficient 
resolution  to  give  distinct  spectral  amplitudes  for  phases 
separated  by  less  than  several  seconds,  even  when  the  phases 
are  distinguished  by  separate  t  .  This  is  illustrated  in 
Section  IV  where  several  synthetic  seismograms  are  analyzed. 
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Figure  6.  Forty-eight  narrow  band  filters  with  center  fre¬ 
quencies  0.2£fjcsl.6  Hz  are  plotted  versus  frequency. 
At  the  bottom  the  filters  for  fk  =  0.2,  0.6,  1.0 
and  1.4  Hz  are  plotted  separately.  Note  that  the 
filter  shape  is  constant  for  fk>0.35,  where  Q  is 
proportional  to  fk. 
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The  variable  frequency  magnitude  (VFM)  is  denoted  by 
For  a  specified  envelope  peak  this  is  defined  by 
(Savino  and  Archambeau,  1974;  Savino  et  al,  1975) 

^(fk)  =  log  (Ak  •  fk)  +  Bp  (A)  ,  (10) 

where  Bp(A)  is  the  usual  Gutenberg-Richtgr.. distance- correction 
for  zero-to-peak  amplitude.  That  is,  B  (A)  is  0.3  uni's  more 

i. 

than  the  peak-to-peak  value  used  in  (1).  This  formula  is^the 
frequency  domain'  analog  of  the  standard  time  domain  formula 
for  m^,  equation  (1) .  Specification  of  the  particular  envelope 
peak  is  based  on  the  t  associated  with  that  peak.  That  is, 
we  choose  the  (narrow)  portion  of  the  wavetrain  from  which  the 
mb(fk)  is  determined. 

3.3  DEFINITION  OF  mb 

/v 

The  motivation  for  the  definition  of  mfa  is  the  observa¬ 
tion  that  the  VFM,  provides  an  estimate  of  the  spectral 

energy  in  a  narrow  band  in  time  and  frequency.  The  idea  is  to 
compute  some  average  of  this  quantity  in  the  region  of  maximum 
spectral  energy  near  the  P  wave  arrival  time  which  is  also  the 
region  sampled  by  conventional  time  domain  m^  measurements. 

/\ 

We  illustrate  the  m,  determination  with  the  HNME  record- 

b 

ings  from  Figure  3.  These  records  were  processed  by  MARS  and 
the  results  are  plotted  in  Figure  S.  These  plots  were  obtained 
in  the  following  way: 

•Each  seismogram  was  Fourier  transformed  and  the 
instrument  response  was  removed  from  the  spectrum. 

•Each  seismogram  was  filtered  by  forth-eight  constant 
width  narrow  band  filters  with  center  frequencies  (fk) 
from  0.2  -  1.6  Hertz  and  Q  specified  by  Equation  (5). 
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•Envelope  peak  amplitudes,  A^(fk),  with  group  arrival 

times,  t  {f.  ) ,  within  an  eight  second  window  about 
g  k 

the  P  wave  arrival  time  were  selected.  In  many  cases 
there  are  two  or  more  such  peaks. 

•  The  largest  peak  for  0.35sfj^l.6  was  identified. 

Using  this  peak  as  a  starting  point,  one  peak  is 
then  selected  for  each  f^,  with  the  selection  made 
to  maximize  the  smoothness  of  t  versus  f^  and  A^ 
versus  f^. 

In  Figure  8  we  show  log  (A^  •  f^)  and  t  versus  f^  for 
each  event.  The  zero  crossing  between  the  first  peak  and  first 
trough  is  defined  to  be  at  t  =  0.  The  amplitude  quantity 
plotted  in  Figure  8  is  essentially  the  VFM,  see  Equation  (3) . 

The  log  (Ak  •  f  )  is  a  velocity-like  quantity.  Since 
the  instrument  response  has  been  removed  and  the  geometric 
attenuation  is  nearly  independent  of  frequency,  the  only  im¬ 
portant  filter  applied  to  the  data  should  be  the  Q  operator; 

that  is,  exp  {  -  irf  t*  }  where  t*  is  the  travel  time  divided 

c 

by  the  effective  Q  for  the  path.  Thus,  log  (Ak  •  f^)  falls 
off  at  high  and  low  frequency  and  has  the  peaked  character 
seen. 

The  data  plotted  in  Figure  8  are  used  to  define  m^. 

There  are  several  possibilities  for  how  this  might  be  done. 

A 

We  define  according  to 

rru  =  log  A  +  B  (A)  +  S  ,  (7) 

b  p  c 


where  B^(A)  is  the  Gutenberg-Richter  distance  correction  and 


Sc  may  be  included  as  some  empirically  determined  station 


correction.  Possible  definitions  for  A  include: 
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1.  A  is  •  fj,  at  a  single  fixed  frequency; 


2.  A  is  the  maximum  value  of  A^  ■  f^  in  a  particular 
f^  band; 


3.  A  is  the  average  A^  •  f^  over  a  fixed  band; 


4. 


A  is  the  average  value  of  A^ 
defined  by  requiring  that  t 

y 

specified  limits. 


•  fk  in  an  £k 
remain  within 


band 


These  alternatives  have  been  evaluated  by  testing  their 
application  to  the  HNME  recordings  of  Pahute  Mesa  explosions 
described  in  Section  II.  These  data  lead  to  the  conclusion 
that  2  and  4  are  the  best  definitions.  They  seem  to  give 
equally  good  results.  We  are  wary  of  recommending  any  def¬ 
inition  based  on  a  single  frequency  value  and  are  biased 
toward  the  definition  requiring  some  averaging  over  a  frequency 
band.  Thus,  we  prefer  the  tg  band  —  averaged  A,  definition  4. 

In  Figure  8  the  maximum  value  of  log  (A^  •  f^)  for 

f^  -  0.35  is  indicated  with  an  asterisk.  From  this  point  our 

algorithm  searches  backwards  and  forwards  in  f^  to  define  the 

limits  of  the  region  within  which  the  tg  are  within  one  second 

of  the  t  at  the  peak.  The  limits  of  this  band  are  indicated 

q 

by  small  vertical  bars  on  the  plots.  Then  A  is  defined  by 


* 


A 


max 


min 


(Ak  •  fk)df 


(8) 


where  Af^  =  max  f^  -  min  f^,  with  the  latter  two  values  being 
those  indicated  on  the  plots.  The  integral  is  evaluated  by 
quadratures.  Then  m^  is  computed  according  to  (7). 


<# 
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The  scheme  we  have  outlined  is  easily  automated  as  a 
post-processor  of  the  MARS  output  table  of  A^  -  t  pairs  as  a 
function  of  f.  .  As  input  we  select  a  time  window  specifying 
the  portion  of  the  signal  where  we  want  m^  to  be  determined. 
Envelope  peaks  with  t  outside  this  window  are  discarded. 

For  the  seismograms  analyzed  in  this  report  this  window  ex¬ 
tended  from  three  seconds  before  to  five  seconds  after  the 
first  zero  crossing  (the  t  =  0  on  the  plots  of  Figure  8) . 

We  also  specify  a  window  for  the  f^  and  require  the  limits 
for  the  integration  (8)  to  be  within  that  window.  For  this 
analysis  no  upper  limit  for  f^  was  necessary.  The  lower  limit 
was  usually  set  at  f^  =  0.35  Hz.  This  limit  was  set  higher  in 
cases  where  the  low  frequency  log  (A^  •  f^)  seemed  to  be  noise 
contaminated  and  we  will  point  this  out  at  the  appropriate 
place. 

In  Figure  9  we  again  plot  the  HNME  seismograms  from 
Figure  4.  A  time  band  is  indicated  on  each  seismogram.  The 
tg  for  all  Ak  used  in  the  computation  of  mfa  fall  within  this 
band.  The  tg  associated  with  the  peak  log  (Afc  •  f^)  is  marked 
with  an  asterisk. 

For  each  t  in  the  indicated  band  the  time  domain  filters 

g 

have  the  form  shown  in  Figure  7.  The  portion  of  the  signal  that 

A 

influences  the  cannot  be  specified  very  precisely,  but 
Figures  7  and  9  give  a  good  indication  of  the  time  window  that 
might  be  important.  In  general,  this  window  is  some  5  to  8 
seconds  long. 

Another  way  to  better  understand  the  sensitivity  of 
is  to  apply  the  algorithm  to  synthetic  seismograms.  This  is 
done  in  Section  IV  for  synthetics  that  are  nearly  identical  to 
the  observed  HNME  recordings  of  MAST,  CAMEMBERT,  and  FONTINA. 


The  HNME  recordings  from  Figure  4  are  shown  with 
a  scale^ indicating  the  limits  of  the  tg  band  use 
in  the  m,  calculation. 


Figure  9 
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m.  FOR  HNME  SEISMOGRAMS 
b 

Using  the  definition  given  in  the  previous  section,  we 
compute  the  m.  from  the  data  plotted  in  Figure  8 .  The  results 
are  summarized  in  Table  5.  We  give  the  m^  and  the  interval  of 
integration  for  (8) .  This  interval  is  also  indicated  by  ver- 

A 

tical  bars  on  the  event  plots  in  Figure  8.  The  m^  is  computed 
from  (7)  with  B  (A)  =  3.53.  This  is  the  same  formula  used  to 
compute  the  time  domain  m^  listed  in  Table  4  and  these  m^  are 
listed.  Also  listed  is  the  frequency  (marked  with  an  asterisk 
in  Figure  8)  for  each  event  and  the  in^  obtained  from  (7)  with 
A  being  the  maximum  log  (A^  •  f^) . 

If  the  shape  of  log  (A^  •  f^)  were  the  same  for  all 
events,  the  difference  between  the  m^  and  mfa  listed  in  Table  5 
would  be  constant.  In  fact,  the  mean  residual  is  0.09  units 
with  a  standard  deviation  of  0.03.  Thus  the  difference  is 


fairly  constant,  though  the  interval  of  integration  varies 
from  event- to-event. 
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•  TABLE  5.  FOR  HNME  SEISMOGRAMS 


fk 

Sb 

min  f 

A 

Table 

Event 

at  peak 

(fk  at  peak) 

max 

“b 

"b 

Stilton 

0.78 

5.89 

0.40-1.05 

5.84 

5.86 

Pool 

0.65 

6.39 

0.45-0.95 

6.28 

6.29 

Estuary 

0.55 

6.34 

0.35-0.88 

6.25 

6.07 

Tybo 

0.60 

6.41 

0.35-0.90 

6.33 

6.26 

Mast 

0.78 

6.39 

0.35-1.13 

6.24 

6.25 

Cheshire 

0.45 

6.31 

0.45-0.83 

6.23 

6.02 

Camembert 

0.55 

6.69 

0.35-0.80 

6.60 

6.37 

Muenster 

0.40 

6.70 

0.35-0.75 

6.63 

6.49 

Colby 

0.50 

6.88 

0.35-0.75 

6.78 

6.50 

Kasseri 

0.50 

6.90 

0.35-0.80 

6.83 

6.59 

Fontina 

0.40 

6.91 

0.35-0.75 

6.79 

6.43 
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IV.  mb  FOR  SYNTHETIC  SEISMOGRAMS 

Considerable  insight  into  the  behavior  of  the  MARS 
filters  and  their  use  to  define  in^  can  be  gained  by  studying 
synthetic  seismograms .  In  this  section  we  first  construct 
synthetic  seismograms  at  the  LRSM  station  HNME  for  three 
Pahute  Mesa  events,  MAST,  CAMEMBERT  and  FONTINA.  We  then 

A 

process  these  synthetic  records  with  the  m^  algorithm. 

The  synthetic  seismograms  are  constructed  in  the 
following  way: 

•The  computational  method  is  that  described  by  Bache 
and  Harkrider  (1976). 

•The  explosion  source  function  is  a  reduced  displacement 
potential  ('i'(t))  computed  with  the  Mueller/Murphy  the¬ 
ory  (Mueller  and  Murphy,  1971) .  The  amplitude  of  the 
transformed  V (t)  is  shown  in  Figure  10  for  the  three 
events.  The  depth  at  which  the  source  function  was 
computed  is  indicated  on  the  figure. 

•The  crustal  models  for  the  source  region  are  tabulated 
in  Table  6.  Information  used  to  construct  these  models 
included  well-log  data,  average  overburden  velocity 
data  and  Basin  and  Range  crustal  structure  information 
from  refraction  profiles  (Hill  and  Pakiser,  1967) . 

Some  adjustments  to  the  velocity  of  layers  above  the 
source  were  also  made  to  change  the  P-pP  delay  time  to 
improve  the  agreement  of  synthetic  and  observed  records. 

•For  the  crust  in  the  receiver  region  we  used  the  model 
tabulated  in  Table  7.  The  upper  mantle  model  is  HNME 
(Helmberger  and  Wiggins,  1971). 
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Figure  10.  The  source  functions  for  three  Pahute  Mesa  explo¬ 
sions  are  plotted  with  the  amplitude  axis  scaled 
to  0.02  KT  while  the  frequency  axis  is  scaled  to 
600  KT. 
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TABLE  6.  SOURCE  REGION  CRUSTAL  STRUCTURE  FOR 
SYNTHETIC  SEISMOGRAM  CALCULATIONS 


Depth 

(km) 

Thickness 

(km) 

ct 

(km/sec) 

6 

(km/sec) 

P 

(gm/cm 

MAST 

0.32 

0.32 

2.38 

1.30 

1.66 

0.67 

0.35 

3.24 

1.90 

2.23 

1.10 

0.43 

3.90 

2.30 

2.10 

1.58 

0.48 

4.80 

2.60 

2.65 

1.72 

0.14 

4.30 

2.40 

2.55 

1.77 

0.05 

3.85 

2.20 

2.50 

2.10 

0.32 

4.40 

2.50 

2.58 

6.00 

3.90 

4.70 

2.60 

2.60 

12.00 

6.00 

5.40 

2.70 

2.70 

20.00 

8.00 

6.00 

3.50 

2.80 

CAMEMBERT 


0.34 

0.34 

2.92 

1.69 

2.00 

0.42 

0.08 

4.43 

2.56 

2.12 

1.50 

1.08 

3.50 

2.02 

2.10 

2.10 

0.60 

4.30 

2.40 

2.60 

Layers  8- 

■10 

of  MAST  structure 

FONTINA 

Layers  1-2 

of 

CAMEMBERT 

structure 

0.76 

0.34 

2.69 

1.55 

1.93 

0.91 

0.15 

2.88 

1.66 

2.02 

1.50 

0.59 

3.31 

1.91 

2.20 

Layers  4-7 

of 

CAMEMBERT 

structure 

TABLE  7.  RECEIVER  REGION  CRUSTAL  MODEL 


Depth 

(km) 

Thickness 

(km) 

a 

(km/sec) 

8 

(km/sec) 

P 

(km/sec) 

1.7 

1.7 

4.0 

2.31 

2.3 

3.0 

1.3 

5.1 

2.94 

2.5 

20.0 

17.0 

6.0 

3.50 

2.8 
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•  Anelastic  attenuation  was  included  by  the  operator 
(Strick,  1970) . 


with  t*  =  0.8 

•The  response  of  the  KS  36000  seismometer  was  included 
by  a  frequency  domain  operator.  The  amplitude  and 
phase  response  are  tabulated  in  Table  8.  Woolson  (1978) 
states  that  there  is  some  ambiguity  in  the  definition 
of  the  response  of  this  instrument,  in  part  due  to  the 
fact  that  it  is  not  a  minimum  phase  filter.  We  point 
this  out  as  a  possible  source  or  error. 

We  have  accumulated  much  experience  with  comparing 
synthetic  seismograms  computed  this  way  with  observed  seismo¬ 
grams  (e.g.,  Bache  et  aJL.  ,  1975;  Bache  et  al .  ,  1976;  Bache, 
1977) .  In  most  of  our  previous  work  the  source  was  a  spheri¬ 
cally  symmetric  point  source  in  a  layered  elastic  medium. 

The  important  phases  are  then  P  and  pP  which  is  the  P  wave 
times  the  elastic  free  surface  reflection  coefficient.  There 
is,  however,  a  serious  question  about  the  accuracy  of  this 
representation  of  pP.  Comparison  of  synthetic  and  observed 
seismograms,  especially  for  a  wide  range  of  yield,  suggest 
that  the  actual  effect  of  pP  is  less  than  predicted  when  an 
elastic  pP  is  assumed. 

The  explosion  source  is  clearly  much  more  complicated 
than  a  spherically  symmetric  point  source.  Large  amounts  of 
surface  spallation  are  known  to  occur.  Using  data  from  near 
source  gauges,  there  have  been  several  attempts  to  estimate 
the  extent  of  the  spalled  region  (Eisler  and  Chilton,  1964; 


c 


TABLE  S.  INSTRUMENT  RESPONSE  FOR  THE  KS  36000  SEISMOMETER 

(Woolson,  Private  Communication) 

Frequency  Normalized  Amplitude  Phase  (radians) 


0.01 

0.59  x  10"5 

-2.175 

0.02 

0.65  x  10"4 

-2.178 

0.05 

0.41  x  10"3 

-2.194 

0.10 

0.265  x  10-2 

-2.279 

0.20 

0.0165 

-2.458 

0.40 

0.  104 

-2.777 

0.50 

0.204 

-2.984 

0.57 

0.287 

-3.157 

0.67 

0.43 

-3.393 

0.80 

0.66 

-3.707 

1.00 

1.00 

-4.156 

1.  11 

1.09 

-4.325 

1.25 

1.  45 

-4.558 

1.50 

1.80 

-4.983 

1.67 

2.00 

-5.238 

2.00 

2.25 

-5.669 

2.22 

2.35 

-5.927 

2.50 

2.38 

-6.225 

3.00 

2.30 

-6.699 

3.33 

2.11 

-6.976 

4.00 

1.63 

-7.396 

5.00 

1.09 

-7.766 
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Viecelli,  1973;  Sobel,  1978).  We  shall  see  that  the  spall 
impact  phase  resulting  from  the  estimates  made  by  these  authors 
is  large  enough  to  be  clearly  observable  on  the  teleseismic 
records . 

The  synthetic  seismogram  calculations  were  done  with  a 
composite  source  including  the  reduced  displacement  potential 
for  the  explosion  and  a  downward  impulse  applied  at  the  surface 
to  represent  the  impulse  generated  by  spall  closure.  The 
presence  of  large  amounts  of  spallation  also  implies  that  the 
pp  is  not  the  same  as  it  would  be  if  the  material  were  behaving 
elastically.  We  expect  some  degradation  of  pP;  some  loss  of 
energy  from  this  phase.  Here  we  arbitrarily  reduce  pP  by  some 
constant  (frequency- independent)  factor. 

Fixing  the  path  models  as  we  have  described,  we  adjusted 
the  source  to  achieve  a  good  fit  to  the  observed  seismograms. 
The  free  parameters  that  were  adjusted  are  as  follows: 

1.  The  P  -  pP  lag  time  was  adjusted  within  narrow 
limits  imposed  by  the  explosion  depth,  the  known 
overburden  velocities  and  the  observed  source- 
surface  travel  time.  This  was  done  by  adjusting 
the  velocities  of  the  overburden  layers  and/or 
computing  with  a  source  depth  that  is  slightly 
different  than  the  actual  depth. 

2.  The  upgoing  waves  from  the  source  were  multiplied 
by  a  constant  y  <  1. 

3.  A  spall  impulse  given  by  an  amplitude  and  time  lag 
(with  respect  to  the  explosion)  was  added.  These 
two  parameters  were  adjusted  within  limits  imposed 
by  the  empirical  estimates  of  Viecelli  (1973)  and 
Sobel  (1978) .  The  resulting  teleseismic  phase  is 
denoted  P  . 


Our  final  synthetic  seismograms  are  shown  in  Figure  11. 
The  construction  of  these  records  is  illustrated  by  separate 
plots  of  P  (all  upgoing  source  energy  was  deleted  from  the 
calculation)  and  P  +  pP.  The  P  phase  is  shown  for  only  the 

S 

MAST  event.  Since  the  source  is  an  impulse,  the  P  is  the 

b 

impulse  response  of  the  layered  earth  model  and  this  is  nearly 
the  same  for  all  three  events. 

The  comparison  between  synthetic  and  observed  seismo¬ 
grams  is  shown  in  another  way  in  Figure  12.  The  agreement  is 
remarkably  good,  especially  considering  the  extremely  simple 
model  for  spall  impact  used  in  the  calculation. 

The  choices  for  the  free  parameters  used  for  the 
synthetics  in  Figures  11  and  12  are  summarized  in  Tables  9 
and  10  where  they  are  compared  to  independent  estimates  of 
these  parameters.  First,  we  consider  the  P-pP  lag  time  data 
in  Table  9.  The  observed  values  for  mean  overburden  velocity 
(measured  with  small  amplitude  seismic  waves)  and  source-to- 
surface  travel  time  (measured  at  shot  time)  are  not  consistent. 
Both  are  measures  of  the  velocity  of  very  high  frequency  waves. 
The  most  we  would  like  to  say  about  the  P-pP  lag  time  used  in 
the  calculations  is  that  it  is  not  inconsistent  with  the  near¬ 
source  data  summarized  in  Table  9. 

The  parameters  for  the  spall  impulse  are  summarized  in 
Table  10.  We  see  that  the  amplitude  is  between  the  estimates 
of  Viecilli  (1973)  and  Sobel  (1973),  which  is  a  reasonable 
place  to  be.  The  delay  time  we  found  necessary  to  match  the 
data  is  a  bit  shorter  than  the  estimates  given  by  these 
authors,  but  not  by  too  much. 

An  is  given  for  each  seismogram  in  Figures  11  and 
12  except  the  Pg  record.  This  m^  was  determined  with  the  semi- 
automated  time  domain  procedure  described  in  Section  II.  The 
phase  measured  is  first  trough  to  second  peak  in  each  case. 
Concerning  this  m^,  we  see  the  following: 
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Figure 


12.  Comparison  of  synthetic  (heavy  lines)  and 
observed  seismograms. 


TABLE  9.  PARAMETERS  FOR  THE  P-pP  LAG  TIME  IN  THE 
SYNTHETIC  SEISMOGRAM  CALCULATIONS 


Theoretical  parameters 

MAST 

CAMEMBERT 

FONTINA 

Depth  of  burial  (km) 

1.07 

1.10 

1.18 

P-pP  lag  (sec) 

0.67 

0.63 

0.76 

Mean  overburden  velocity  (km/sec) 

3.10 

3.35 

2.99 

Mean  overburden  density  (gm/cm3 ) 

2.01 

2.07 

2.04 

Source-surface  travel  time  (sec) 

0.34 

0.33 

0.39 

Observed  Data  (Tables  1-3) 


Actual  depth  (km) 

0.91 

1.31 

1.22 

Mean  overburden  velocity  (km/sec) 

3.88 

2.90 

2.86 

Mean  overburden  density  (gm/cm3) 

2.24 

2.10 

2.00 

Source-surface  travel  time  (sec) 

0.36 

0.45 

- 

Depth/mean  overburden  velocity  (sec) 

0.23 

0.45 

0.43 
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TABLE  10.  PARAMETERS  FOR  THE  Pe  PHASE  IN  THE 

s 

SYNTHETIC  SEISMOGRAM  CALCULATIONS 


MAST 

CAMEMBERT 

FONTINA 

m 

Spall  Impulse 
(Dyne- sec/km) 

10  x  1014W 

14 

14  x  10  W 

9  x  1014t 

P-Ps  lag  time 
(sec) 

1.25 

1.92 

1.86 

Y* 

0.60 

0.50 

0.50 

Estimated  of  spall  impulse  from  near-field  data 

Viecelli  (1973):  4.6  x  1014W 

Sobel  (1978):  21-25  x  1014W 


Estimated  of  P-P  delay  from  near-field  data: 

s 

2. 0-2. 5  seconds 


♦Factor  multiplying  the  upgoing  waves  from  the  source. 
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•  The  pP  phase  enhances  the  m^  compared  to  that  for  P 
alone.  This  enhancement  is  less  than  would  occur  with 
these  P-pP  lag  times  if  we  had  not  reduced  the  size  of 
the  pP  phase  by  about  half. 

•The  Pg  phase  has  almost  no  effect  on  the  time  domain 
mb‘ 

•While  this  is  not  the  purpose  here,  we  note  that  the 
agreement  of  observed  and  synthetic  m^  is  quite  good. 

This  gives  confidence  in  our  source  and  path  models. 

In  this  report  we  do  not  want  to  make  too  much  of  the 
interpretation  of  the  prominent  later  phase  as  being  generated 
by  spall  closure.  There  are  alternative  explanations  for  this 
secondary  phase.  From  the  data  at  one  station  we  cannot  be 
sure  it  is  not  a  multipathing  phenomenon.  Also,  there  are  other 
source  effects  that  could  generate  such  a  phase.  Tectonic 
strain  release  is  one  that  comes  to  mind.  However,  our  pre¬ 
vious  work  (Bache,  1976)  indicates  that  the  tectonic  release 
component  can  only  be  this  large  for  the  most  favorable  source 
orientation-station  azimuth  combinations.  Again,  we  must  look 
at  data  from  other  stations  to  see  if  tectonic  strain  release 
is  a  plausible  explanation. 

Our  objective  here  is  to  generate  synthetic  seismograms 
that  closely  resemble  the  data  so  we  can  better  understand  our 

A 

algorithm.  Using  the  spall  impulse  model  to  generate  a 
phase  with  the  right  time  delay  and  amplitude,  we  have  ac¬ 
complished  that  objective.  More  detailed  discussion  of  the 
physical  nature  of  this  phase  will  be  given  in  a  separate 
report . 
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We  now  show  the  results  of  applying  our  m^  algorithm 
(Section  III)  to  the  synthetic  seismograms  of  Figure  11.  The 
plots  of  log  (A^  •  f^)  and  t  versus  f^  are  shown  in  Figure  13 
for  the  synthetic  seismograms  and  for  the  observations.  The 
algorithm  used  to  obtain  these  plots  is  described  in  Section  3.3, 

The  calculation  of  m^  from  the  log  (A^  •  f^)  and  t  data 

is  described  in  Section  3.3.  We  begin  by  obtaining  the  peak 

value  of  log  (A^  •  f^) .  This  is  indicated  by  an  asterisk  on 

the  plot.  We  then  define  the  upper  and  lower  limits  of  an  f^ 

band  about  the  f^  at  the  peak.  These  limits  are  chosen  such 

that  the  t  is  within  one  second  of  the  t  at  the  peak  log 
g  g  c  3 

(A^  •  f^) .  The  limits  chosen  automatically  by  this  algorithm 

are  indicated  on  the  synthetic  P  +  pP  +  Pg  and  observed  seismo¬ 
grams  in  Figure  13  by  small  vertical  bars.  For  the  P  and  pP 
synthetics  for  each  event  the  f^  limits  were  forced  to  be  the 

same  as  for  the  P  +  pP  +  P  .  This  was  done  so  we  could  see 

s 

the  effect  of  adding  the  individual  phases  without  being 
mislead  by  averaging  over  different  bands. 

From  the  data  plotted  in  Figure  13  we  compute  in^  from 


m^  =  log  A  +  3.48 


(A.  1 ) 


where 


max  f , 


min  f , 


V  dfk 


(A. 2) 


and  =  max  f^  -  min  f^,  the  band  indicated  on  the  plots. 
On  each  of  the  nine  plots  in  Figure  13,  we  give  the  m^  values 
computed  from  the  above  equations. 
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Figure 


In  Figure  14  we  again  show  the  full  (P  +  pP  +  P  ) 

S 

synthetics  with  the  t  time  band  for  the  m^  calculation  indi¬ 
cated  on  each  seismogram.  The  asterisk  denotes  the  t  asso¬ 
ciated  with  the  peak  log  (A^  •  f^) ,  which  is  also  denoted  by 
an  asterisk  in  Figure  13.  The  bars  indicate  the  maximum  and 

minimum  t  within  the  band.  We  see  that  the  identified  ar- 
9 

rival  time  for  the  energy  contributing  to  m^  is  near  the  front 
of  the  record  where  conventional  time  domain  m^  measurements 
are  made.  However,  since  the  width  of  the  filter  is  about 
four  seconds  (defined  as  ±  one  standard  deviation,  see  Figure 
10  in  Section  3.3),  energy  arriving  outside  the  band  indicated 

A 

in  Figure  14  contributes  to  the  m^. 

/v 

The  mb  values  from  Figure  11  and  the  m^  from  Figure  13 
are  summarized  in  Table  11.  From  the  data  in  this  table  and 
the  figures,  we  draw  the  following  conclusions  about  m^: 

•The  agreement  of  observed  and  synthetic  m,  is  remarkable, 

A  ^ 

and  that  for  is  even  better. 

A 

•Comparing  P  and  P  +  pP,  the  m^  and  differences  are 
the  same.  Thus  the  pP  phase  seems  to  have  about  the 
same  effect  on  both  magnitude  measures,  at  least  for 
this  depth  range. 


•  The  P  phase  has  no  effect  on  m.  .  However,  it  enhances 
As  b 

the  by  about  0.10  units. 

•Considering  the  width  of  the  time  domain  filters  and 
the  t  associated  with  the  amplitudes  used  to  compute 

A  /\ 

m^  (Figure  14),  the  effect  of  the  Pg  phase  on  is 
to  be  expected. 
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Figure  14.  The  P  +  pP  +  Pg  synthetic  seismograms  from 
Figure  11  are  shown  with  marks  indicating 
the  t  time  band  for  the  m,  determination. 


TABLE  11. 


SUMMARY  OF  It^ 


AND  m. 


DATA 


A 


mb 

mb 

Seismogram 

(Figure  A. 2) 

(Figure  . 

MAST 

• 

P 

5.39 

5.95 

P  +  pp 

6.03 

6.06 

P  +  pP  +  P 

s 

6.04 

6.15 

Observed 

6.13 

6.19 

CAMEMBERT 


P 

6.12 

6.34 

P  +  pP 

6.21 

6.42 

P  +  pP  +  P 

s 

6.21 

6.54 

Observed 

FONTINA 

6.25 

6.55 

P 

6.30 

6.51 

P  +  pP 

6.38 

6.62 

P  +  pP  +  P 

s 

6.38 

6.72 

<0 

Observed 

6.31 

6.74 
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•From  the  plots  in  Figure  13  we  can  infer  the  signifi¬ 
cance  of  the  spectral  holes  in  the  log  (Ak  •  fk)  plots 
for  the  observed  seismograms.  They  are  almost  cer¬ 
tainly  due  to  the  presence  of  Ps«  The  resolution  of 
the  narrow  band  filters  is  not  nearly  fine  enough  to 
see  effects  of  the  P  -  pP  interference. 
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V.  THREE-DIMENSIONAL  FINITE  DIFFERENCE  EARTHQUAKE  MODELING 

5.1  INTRODUCTION 

In  this  section,  we  summarize  and  discuss  the  results 
of  our  earthquake  modeling  with  the  ILLIAC,  using  TRES,  a 
three-dimensional  finite  difference  code  (Cherry,  1977) .  We 
modeled  faulting  as  a  propagating  stress  relaxation  due  to 
shear  failure  on  a  planar  surface.  Ultimately,  we  would  like 
to  specify  the  relevant  physical  properties  of  the  medium 
and  its  initial  conditions,  and  allow  a  mathematical  model 
of  failure  to  determine  the  subsequent  evolution  of  the  fault 
plane.  However,  at  this  stage,  we  have  not  addressed  the 
physical  mechanism  of  failure.  Instead,  we  prescribe  the 
propagation  of  the  fault  surface,  and  boundary  conditions 
on  the  fault  surface  are  governed  by  a  simple  Coulomb 
friction  law. 

VJhile  the  TRES  algorithm  is  quite  flexible  in  the 
allowed  specification  of  geometry  and  material  behavior, 
the  version  currently  operational  on  the  ILLIAC  is  somewhat 
restricted.  The  faulting  must  nucleate  from  a  point  and 
propagate  with  circular  symmetry  until  reaching  the  edges 
of  a  rectangular  fault  plane.  There  are  no  material  bound¬ 
aries  other  than  the  fault  plane;  that  is,  the  calculations 
are  done  in  a  whole  space.  The  material  behavior  is  linearly 
elastic  except  in  the  vicinity  of  the  fault  plane  where 
plastic  yielding  is  permitted. 

We  performed  two  three-dimensional  finite  difference 
calculations  for  this  study.  Details  of  the  fault  model  are 
given  by  Day,  et  al.  (1978).  The  two  calculations  differed 
only  in  the  yield  strength  Y  assigned  the  material.  Both 
calculations  were  for  a  square  fault  plane  in  a  uniform 
whole  space,  with  rupture  initiated  at  the  center  of  the 
square  fault  (Figure  15) .  The  following  parameters  were 
employed  for  both  calculations: 
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P  wave  velocity 
S  wave  velocity 
Density 

Rupture  velocity 
Tectonic  shear  stress 
Frictional  stress 
Fault  dimensions 


a  =  5.93  km/sec 
3  =  3.42  km/sec 
p  =  2.74  gm/cm 
VR  =  3.08  km/sec 
a-p  =  1  kbar 
Cf  =  .82  k  bars 
2ax2a  =  3  km  x  3  km 


For  the  first  finite  difference  calculation,  the 
"elastic"  model,  the  yield  strength  Y  was  set  to  infinity, 
so  the  constitutive  model  was  linearly  elastic.  For  the 
second  finite  difference  calculation,  the  "plastic"  model, 

Y  was  set  to  /3aT,  so  that  the  fault  zone  was  initially 
stressed  to  the  failure  surface.  With  this  choice  of  Y, 
plastic  flow  ensues  when  the  second  deviatoric  stress 
invariant  J2  increases  above  its  initial  value.  This 
dissipates  any  dynamic  shear  stress  concentration  ahead  of 
the  crack  tip.  Plastic  yield  was  permitted  only  within 
0.2  km  of  the  fault  plane;  elsewhere  linear  elasticity 
was  employed. 

For  both  the  elastic  and  plastic  fault  problems, 
the  medium  was  represented  by  cubic  finite  difference  zones 
0.1  km  on  a  side.  The  numerical  grid  was  large  enough  that 
no  reflection  from  the  exterior  grid  boundary  returned  to 
the  fault  zone  during  the  calculation. 

These  initial  calculations  permit  us  to  examine  the 
accuracy  of  the  three-dimensional  numerical  method,  investi¬ 
gate  the  near-  and  far-field  signal  froma  a  simple  propagating 
stress-relaxation  model,  and  clarify  the  physical  interpreta¬ 
tion  of  the  parameters  of  the  Archambeau  spherical  source 
model.  The  inclusion  of  plastic  yielding  in  the  second  cal¬ 
culation  is  a  sisnificant  step  toward  incorporating  realistic 
rock  mechanics  into  the  fault  model,  and  we  examine  its  effect 
on  the  near-  and  far-field  signals. 
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5.2 


SUMMARY  OF  RESULTS 


The  main  results  of  these  two  calculations  are  sum¬ 
marized  here;  for  details,  we  refer  to  Day,  et  al . ,  (1978). 

5.2.1  Comparison  of  an  Analytic  Solution 

Kostrov  (1964)  obtained  an  analytical  solution  for 
the  problem  of  a  circular  crack  which  nucleates  at  a  point 
in  a  homogeneous,  unbounded  elastic  medium  and  expands  at 
a  constant  rupture  velocity  without  stopping.  This  problem 
corresponds  exactly  to  the  conditions  of  the  elastic  model 
considered  in  this  study,  until  time  a/VR,  where  a  is  the 
half-length  of  the  square  fault.  We  can  compare  the  initial 
fault  slip,  obtained  numerically,  to  Kostrov's  analytic 
solution,  although  once  the  rupture  front  reaches  the  edge 
of  the  square  fault  plane  and  stops  growing,  we  expect  the 
numerical  solution  to  begin  to  deviate  significantly  f.  m 
Kostrov's  solution. 

Figure  16  shows  the  slip  obtained  at  several  points 
in  the  fault  plane.  The  dashed  curves  are  the  finite  dif¬ 
ference  solution  and  the  solid  curves  are  the  analytic 
solution.  The  vertical  bars  indicate  the  arrival  times  of 
edge  effects  due  to  stopping  of  the  rupture  at  its  outer 
boundary.  The  two  solutions  display  the  anticipated  agree¬ 
ment  at  each  point  prior  to  the  arrival  of  the  edge  effects. 
The  small  deviation  of  the  numerical  solution  from  the 
analytic  solution  at  early  time  results,  at  least  partially, 
from  imprecise  weighting  of  the  stress  drop  to  account  for 
the  fractional  rupture  of  a  finite  difference  zone  by  the 
circular  rupture  front.  Archuleta  and  Frazier  (1978) 
achieved  somewhat  better  agreement  with  Kostrov's  solution 
by  incorporating  fractional  rupture  into  their  finite 
element  scheme. 


TIME  (seconds) 


Figure  16.  Relative  displacement  on  the  fault  for  the  elastic 
case.  The  dashed  curves  are  Kostrov's  analytic 
solution;  the  solid  curves  are  the  finite  difference 
results.  x,y  coordinates  in  kilometers  are  given  in 
parenthesis.  Vertical  lines  indicate  the  arrival 
times  of  edge  effects  due  to  fault  finiteness. 


5.2.2  Stress  History  Near  the  Fault 

Figure  17  depicts  the  stress  histories  near  the  fault 
plane  for  the  two  models  (the  purely  elastic  model  and  the 
plastic  model) .  The  shear  stress  in  the  direction  of  pre¬ 
stress  (OyZ)  is  plotted  versus  time  for  seven  points  along 
the  fault  plane  diagonal,  at  increasing  distance  from  the 
hypocenter.  The  stresses  shown  are  actually  evaluated  at 
the  finite  difference  zone  centers  adjacent  to  the  fault, 
which  are  0.05  km  from  the  fault  plane. 

First,  consider  the  stress  for  the  elastic  problem 
shown  in  Figure  17.  Initially,  the  stress  at  a  given  point 
is  at  the  prestress  level  (1  kbar) .  Prior  to  the  rupture 
front  arrival  at  a  given  location,  the  stress  increases 
above  the  prestress  level.  This  stress  concentration  ahead 
of  the  rupture  front  is  a  general  characteristic  of  elasto- 
dynamic  cracks  with  subsonic  rupture  propagation.  We  note 
the  amplification  of  this  stress  concentration  with  increas¬ 
ing  distance  in  the  direction  of  rupture.  At  the  rupture 
arrival  time,  the  stress  drops  abruptly,  over-shooting  and 
then  settling  at  the  prescribed  frictional  stress  of  820 
bars.  The  over-shoot  results  from  the  fact  that  the  observa¬ 
tion  points  are  slightly  removed  from  the  fault  plane  itself. 
The  same  phenomenon  is  present  in  Richards  (1976)  analysis 
of  the  self-similar  expanding  crack.  The  stress  then  remains 
at  the  frictional  stress  level  until  the  nearby  part  of  the 
fault  plane  heals.  Then  the  stress  relaxes  to  a  value  less 
than  that  of  kinetic  friction.  The  fault  over-shoots  the 
static  equilibrium  value  of  slip.  The  healing  wave  shows 
up  clearly  in  Figure  17;  it  propagates  toward  the  hypocenter 
from  the  periphery  of  the  fault.  The  last  phase  evident  in 
the  figure  propagates  outward  from  the  hypocenter,  and  this 
phase  corresponds  to  the  shear  wave  associated  with  the 
final  arrest  of  slip  at  the  center  of  the  fault. 

Now  consider  the  stresses  for  the  plastic  problem  in 
Figure  17.  Until  healing  occurs  at  a  given  point,  the  stress 
history  is  unchanged  from  the  elastic  case  except  that  the 
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stress  concentration  ahead  of  the  rupture  front  has  been 
eliminated.  After  healing,  the  residual  stress  is  nearly 
the  same  as  in  the  elastic  case. 


Figure  18  displays  the  stress  <jy2  near  the  fault  plane 
as  a  function  of  azimuth  <J>.  Again,  the  stress  is  zone  cen¬ 
tered,  so  it  is  actually  evaluated  at  points  0.05  km  dis¬ 
placed  from  the  fault  plane.  The  five  points  displayed  are 
at  approximately  the  same  distance,  1.15  km,  from  the  center 
of  the  fault.  The  stress  concentration  proceding  rupture  is 
nearly  zero  (actually  slightly  negative)  at  <f>  =  0°  ,  and  in¬ 
creases  smoothly  to  a  maximum  at  $  =  90°  .  This  pattern  can 
be  compared  to  Figure  8  of  Richards  (1976).  The  stress 
histories  are  very  similar,  although  his  results  are  for  an 
elliptical  fault  in  which  rupture  velocity  varies  with 
azimuth  from  0.92  8  to  0,  whereas  our  numerical  solution  is 
for  circular  rupture  propagation  at  rupture  velocity  0.9  0. 

Figure  19  shows  the  stress  component  ayy  along  the 
fault  plane.  This  component  of  stress  is  not  relieved  by 
plastic  flow,  and  the  concentration  of  ayy  ahead  of  the 
rupture  is  essentially  identical  in  both  the  plastic  and 
the  elastic  cases. 

5.2.3  Velocity  History  on  the  Fault 

Figure  20  shows  the  slip  velocity  obtained  on  the 
fault  plane  at  increasing  distances  from  the  hypocenter 
along  a  radial  line.  The  solid  curves  are  for  the  elastic 
case,  the  dashed  curves  are  for  the  plastic  case. 

It  is  evident  from  Figure  20  that  the  initial  velocity 
is  strongly  peaked  and  the  peak  value  increases  with  hypo- 
central  distance.  We  can  understand  this  characteristic  of 
the  velocity  curves  by  means  of  Kostrov’s  analytic  solution. 
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Figure  18.  Time  history  of  a yZ  adjacent  to  the  fault  plane 

for  several  azimuthal  angles  </>.  Each  observation 
point  is  at  a  distance  0.05  km  from  the  fault 
plane,  and  a  distance  1.15  km  from  the  center  of 
the  fault. 


which  predicts  that  peak  slip  velocity  should  increase  as 
the  square  root  of  the  hypocentral  distance.  This  prediction 
is  in  good  agreement  with  the  numerical  results  in  Figure  20. 

It  is  evident  from  Figure  20  that  the  initial  part  of 
the  particle  velocity  is  unaffected  by  the  introduction  of 
plasticity,  within  the  resolution  of  the  finite  difference 
calculation.  A  large  velocity  peak  occurs  at  the  crack  tip, 
even  in  the  absence  of  the  strong  stress  concentration  asso¬ 
ciated  with  the  purely  elastic  problem.  The  plastic  and 
elastic  solutions  are  indistinguishable  until  the  arrival 
of  the  stopping  phase. 

After  the  arrival  at  a  given  point  of  edge  effects 
due  to  stopping  of  the  rupture  front,  the  fault  plane  veloc¬ 
ities  are  substantially  modified  by  the  plasticity.  As 
Figure  20  indicates,  yielding  at  the  crack  tip  smooths  the 
stopping  phase,  robbing  the  slip  function  of  high  frequencies 
and  increasing  the  long-period  content  of  the  slip  function. 
The  average  static  slip  on  the  fault  plant  is  increased  by 
about  11  percent  when  yield  is  permitted. 

5.2.4  Radiated  Fields 

The  average  slip  for  the  elastic  fault  was  79  cm, 
which  is  0.61  times  the  value  of  slip  at  the  center  of  the 
fault.  The  seismic  moment  obtained  from  this  average  slip 
was  2.28  x  10^  dyne-cm.  This  value  is  14  percent  greater 
than  the  prediction  obtained  by  combining  the  static  circular 
crack  formula  (Keilis-Borok,  1959)  with  the  expression  for 
seismic  moment,  MQ  -  yAs.  This  over-shoot  of  the  static 
solution  has  been  observed  in  previous  dynamic  modeling  — 
notably  the  work  of  Madariaga  (1976) . 

Figures  21  and  22  present  normalised  far-field  P  and  S 
wave  displacement  spectra  and  time  histories  for  the  elastic 
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case.  Results  are  shown  at  10°  intervals  in  0;  Figure  21 
is  for  4>  =  90° ,  and  Figure  22  is  for  <p  =  45°.  Solid  lines 
are  the  P  wave  displacements,  dashed  lines  are  the  S  wave 
displacements.  The  curves  are  normalized  by  the  zero  fre- 
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quency  spectral  amplitude,  derived  from  the  average  slip  s. 

?  o  'V  —1 

For  P  waves,  the  normalization  is  £>  /a  As  ( 4tt  a  r)  Rp. 

For  S  waves,  we  can  interpret  the  curves  either  as  the  0 

^  -1 

component  of  displacement  normalized  by  As  ( 4 tt  6  r)  RsQ  , 

or  as  the  <j>  component  of  displacement  normalized  by 
%  _1 

As  (4ir  6  r)  R  ,  .  R  ,  R  a  and  R_,  are  double  couple  radia 

s<J>  p  S  0  S<p 

tion  patterns  (see  Day,  et  al.  1978)  ,  A  is  the  fault  area, 
and  r  the  hypocentral  distance.  The  travel  times  from  hypo 
center  to  receiver  have  been  removed  from  the  P  and  S  time 
histories . 


Comparing  results  at  <j>  =  90°  with  those  at  <p  =  45°, 
we  note  that  pulse  width  and  corner  frequency  have  practi¬ 
cally  no  dependence  on  4>;  at  higher  frequency  there  is 
some  difference  in  spectral  and  time  domain  detail  between 
the  2  azimuths.  (Results  at  <p  =  0°  are  virtually  identical 
to  those  at  <j>  =  90°,  and  are  not  shown.) 

Dependence  of  pulse  width  and  corner  frequency  on  0 
and  on  wave  type  (P  or  S)  is  significant.  Our  observations 
concur  with  those  of  Madariaga  (197  6)  : 


(a)  S  wave  corner  frequencies  are  smaller  than  P 
wave  corner  frequencies,  except  near  0=0°. 

(b)  Pulse  width  and  corner  frequency  are  governed 
by  the  travel  time  difference  between  stopping 
phases  from  the  near  and  far  edges  of  the  fault 
Thus,  pulse  width  increases  with  0,  being 
greatest  for  observers  near  the  plane  of  the 
fault  and  smallest  for  observers  near  the 
fault  normal. 
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(c)  P  and  S  wave  corner  frequencies,  for  9  >30°, 
are  expressed  very  well  by  Madariaga’s 
Equation  (24) ,  replacing  the  circular  fault 
radius  with  the  square  fault  half-width. 

The  average  static  slip  for  the  plastic  fault  problem 
is  88  cm,  which  exceeds  the  average  slip  in  the  elastic  case 
by  11  percent.  This  results  from  the  smoothing  of  the  stop¬ 
ping  caused  by  yielding  at  the  fault  edge.  Actually,  this 
is  precisely  the  percentage  increase  that  would  be  predicted 
by  simple  scaling  of  the  elastic  problem  for  a  fault  dimen¬ 
sion  occupying  the  entire  length  and  width  of  the  plastic 
zone  (the  plastic  zone  extended  0.2  km  beyond  edge  of  the 
fault)  . 


The  seismic  moment  for  the  plastic  problem  is  3.15  x 
10 2 4  dyne-cm,  which  is  38  percent  larger  than  the  moment  for 
the  elastic  case.  Scaling  of  the  elastic  solution  as  sug¬ 
gested  in  the  last  paragraph  would  predict  a  slightly  larger 
increase  in  moment,  42  percent  instead  of  38  percent. 

Figure  23  shows  the  effect  of  plastic  yield  on  the 
far-field  displacements.  Spectra  and  pulses  are  shown  at 
<p  =  90°  for  3  values  of  0 :  S  wave  solutions  are  shown  at 
9=0°  and  9  =  90°,  and  P  wave  solutions  are  shown  at  9  =  45 
Dashed  curves  are  the  elastic  case,  solid  curves  the  plastic 
case.  In  each  case,  the  far-field  solution  is  the  sum  of 
multipolar  terms  up  to  £  =  8. 

The  main  influence  of  plastic  yielding  is  to  smooth 
the  stopping  phases,  with  the  result  that  the  low-frequency 
part  of  the  spectrum  is  enhanced  at  the  expense  of  the  high 
frequency  part  of  the  spectrum.  Consider  for  example,  the 
P  wave  pulse  at  9  =  45°.  Two  stopping  phases,  corresponding 
to  rupture  arrival  at  the  near  and  far  edges  of  the  fault. 
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respectively,  are  apparent  in  the  P  wave  displacement  pulse 
for  the  elastic  problem.  These  appear  as  discontinuities 
in  slop  occurring  at  about  0.46  seconds  and  0.67  seconds. 

The  P  wave  pulse  for  the  plastic  problem  coincides  with  that 
for  the  elastic  problem  until  the  arrival  of  the  first  stop¬ 
ping  phase.  The  displacement  pulse  for  the  plastic  case 
then  reverses  more  gradually,  over-shooting  the  elastic 
case;  the  second  stopping  phase  is  almost  imperceptible 
for  the  plastic  case. 

Clearly,  the  increase  in  seismic  moment  is  the  conse¬ 
quence  of  plastic  strain  induced  at  the  periphery  of  the  slip 
surface.  Relieving  stress  by  permitting  plastic  strain  is 
very  similar  to  relieving  stress  by  permitting  frictional 
sliding,  but  with  the  frictional  stress  approximately  equal 
to  Y//3  (which  in  this  case  equals  the  prestress  aT) .  Out¬ 
side  the  zone  in  which  yielding  was  permitted,  a  static  shear 
stress  concentration  about  23  percent  in  excess  of  the  pre¬ 
stress  developed.  Thus,  if  a  larger  plastic  zone  had  been 
specified,  the  seismic  moment  would  have  been  even  greater, 
as  plastic  strains  extended  outward  to  eradicate  the  stress 
concentration.  On  the  other  hand,  had  a  somewhat  higher 
yield  strength  been  specified,  the  results  of  the  plastic 
problem  would  have  approached  those  of  the  elastic  case. 

If  the  yield  papameter  Y//3  had  exceeded  the  tectonic  stress 
by  1.44  (aT  -  a (that  is,  if  Y  had  been  26  percent  larger), 
no  yielding  would  have  occurred,  and  the  two  solutions  would 
have  been  indistinguishable. 

5.2.5  Comparison  to  Radiation  from  an  Archambeau-Minster  Source 

As  a  generator  of  teleseismic  signals,  the  Archambeau- 
Minster  source  model  possesses  several  of  the  relevant  features 
expected  of  earthquake  sources.  We  compared  synthetic  short- 
period  teleseismic  P  waves  generated  by  a  bilateral  version 
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of  the  Archambeau-Minster  source  model  with  those  generated 
by  the  elastic  fault  model.  The  bilateral  Archambeau  source 
consists  of  two  independent,  adjacent  spherical  sources, 
propagating  in  opposite  directions.  For  the  comparison, 
both  models  were  scaled  to  have  the  same  seismic  moment, 
and  the  cross-sectional  area  of  the  Archambeau-Minster 
spherical  volumes  was  made  equal  to  the  area  of  the  square 
fault  plane  (9  km).  Thus,  the  maximum  source  dimension  was 
somewhat  larger  for  the  spherical  source. 

In  Figure  24  we  compare  five  short  period  teleseisms 
from  the  elastic  finite  difference  calculation  to  those  from 
a  bilateral  version  of  the  Archambeau  model.  The  amplitude 
(corrected  for  period-dependent  instrument  response)  and 
period  data  from  these  synthetic  seismograms  are  summarized 
in  Table  12.  The  Archambeau  model  synthetics  are  somewhat 
longer  period  and  lower  amplitude  than  the  elastic  finite 
difference  source  synthetics.  This  reflects  the  slightly 
higher  corner  frequency  produced  by  the  finite  difference 
source,  as  a  result  of  its  smaller  source  dimension. 

The  physical  meaning  of  the  important  parameter  stress 
drop  is  made  more  clear  after  comparing  the  Archambeau  and 
elastic  finite  difference  models.  In  order  to  scale  the 
Archambeau  source  to  have  the  same  moment  as  the  square 
fault  model,  it  was  necessary  to  make  the  parameter  "stress 
drop"  in  the  Archambeau  model  a  factor  of  3.6  smaller  than 
that  of  the  fault  model.  Thus,  stress  drops  are  under¬ 
estimated  by  about  this  amount  if  they  are  based  on  the 
Archambeau  model. 
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TABLE  12 
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AMPLITUDES  AND  PERIODS  OF  THE  SEISMOGRAMS  OF  FIGURE  24 


P  PHASE 

MAXIMUM  PHASE 

Stations 

Amplitude 

(microns) 

Period 

Amplitude 

(microns) 

Period 

■kp^max 

RES 

EFD  Source* 

2.11 

1.1 

2.59 

1.0 

1.23 

Bi-Model  I 

1.84 

1.3 

1.96 

1.2 

1.07 

DAG 

EFD  Source 

1.38 

1.1 

2.38 

1.2 

1.73 

Bi-Model  I 

1.19 

1.2 

2.12 

1.3 

1.78 

KTG 

EFD  Source 

1.29 

1.1 

2.42 

1.2 

1.88 

Bi-Model  I 

1.11 

1.2 

2.24 

1.3 

2.02 

BOG 

EFD  Source 

0.76 

1.1 

1.75 

1.3 

2.30 

Bi-Model  I 

0.69 

1.2 

1.53 

1.4 

2.21 

pro 

EFD  Source 

0.67 

1.1 

1.96 

1.3 

2.93 

Bi- Model  I 

0.61 

1.2 

1.88 

1.4 

3.08 

*  Elastic  Finite  Difference  Source  Model. 
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5.3 


DISCUSSION 


The  inclusion  of  a  simple  form  of  inelastic  material 
response  was  intended  as  a  step  toward  developing  realistic 
models  of  earthquake  physics.  It  was  found  that  the  initial 
portion  of  the  slip  time-function  is  unaffected  by  the  admis¬ 
sion  of  a  simple  form  of  plasticity.  The  large  velocity  peak 
at  the  crack  tip,  characteristic  of  elastic  crack  problems, 
persisted  in  the  inelastic  case.  The  stopping  phase  was  modi¬ 
fied  somewhat,  however.  Yielding  at  the  edge  of  the  fault 
resulted  in  less  abrupt  stopping,  reducing  the  high-frequency 
content  of  both  the  slip  function  and  the  far-field  displace¬ 
ments,  and  increasing  the  average  slip  by  11  percent.  Accu¬ 
mulation  of  plastic  strain  beyond  the  fault  edge  resulted  in 
a  38  percent  higher  moment  than  in  the  elastic  case. 

These  effects  of  plasticity  depend  upon  our  choice  of 
the  magnitude  of  prestress,  the  yield  strength,  the  frictional 
stress,  and  the  dimensions  of  the  nonlinear  zone.  The  problem 
treated  was  an  extreme  case  in  the  sense  that  the  prestress 
level  everywhere  equaled  the  strength  of  the  medium  (Y//3) . 

With  a  moderate  increase  in  the  yield  strength  (26  percent), 
or  a  similar  decrease  in  the  prestress,  there  would  have  been 
no  yielding,  and  the  solution  would  have  been  identical  to 
that  of  the  elastic  problem.  On  the  other  hand,  we  arbitrarily 
limited  the  extent  of  the  plastic  zone.  Outside  the  plastic 
zone,  a  static  shear  stress  concentration  persisted  which  was 
nearly  as  large  as  that  of  the  elastic  case;  had  the  plastic 
zone  been  larger,  more  plastic  strain  would  have  occurred, 
resulting  in  an  even  larger  seismic  moment. 

The  next  steo  will  be  incorporation  into  the  model  of  a 
fracture  criterion,  so  that  rupture  advance  is  governed  by 
rock  strength  rather  than  being  arbitrarily  prescribed.  There 
are  three  requirements  that  we  shall  impose  on  a  realistic 
fracture  mechanism.  The  first  two  are  that  energy  should  be 
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dissipated  at  the  crack  tip  and  that  stresses  and  velocities 
should  remain  bounded  everywhere.  The  third  requirement  is 
that  the  criterion  of  rupture  can  be  formulated  numerically 
so  as  to  be  nearly  zone-si^a  dependent. 

Ida  (1972)  has  shown  that,  when  the  stress  drop  is 
abrupt,  limiting  the  shear  stress  on  the  fault  plane  is  not 
sufficient  to  prevent  singularities  in  velocity  and  stress 
components.  This  is  supported  by  our  observation  (Sections 
5.2.2  and  5.2.3)  that  plastic  yielding  did  not  reduce  peak 
velocities.  On  the  other  hand,  Ida  (1972)  and  Andrews  (1976) 
have  shown  that  a  slip-weakening  model,  in  which  the  stress 
drops  gradually  as  a  function  of  the  relative  displacement, 
produces  bounded  stresses  and  velocities  and  absorbs  finite 
energy  at  the  rupture  front.  Since  the  rupture  front  is 
smeared  out  in  a  slip-weakening  model,  it  is  likely  to  yield 
a  rupture  advance  which  is  nearly  independent  of  zone  size. 

5 . 4  RESEARCH  PLANS 

•  Extend  our  work  on  faults  with  prescribed  rupture  veloc¬ 
ity  to  include  unilateral  rupture  on  a  long,  narrow  fault.  We 
anticipate  performing  one,  or  perhaps  two,  elastic  calculations 
to  evaluate  the  effect  of  this  geometry  on  the  source  function. 

•  Formulate  a  failure  criterion  in  three-dimensions  and 
develop  an  algorithm  to  be  incorporated  into  the  finite  dif¬ 
ference  code. 

•  Exercise  the  rupture  model  in  the  presence  of  constant 
prestress.  Relationship  of  the  model  parameters  to  rupture 
velocity,  stress  concentrations  and  slip  function  shape  will 
be  examined.  It  is  anticipated  that  rupture  growth  will  have 
to  be  terminated  artificially  in  this  model. 


•  Investigate  the  circumstances  under  which  rupture 
growth  stops  spontaneously.  An  important  mechanism  may  be 
nonuniform  stress  drop.  Several  studies  based  on  kinematic 
earthquake  models  (for  example,  Bache  and  Barker,  1978  and 
Barker,  et  al. ,  1978)  have  indicated  that  variable  stress 
drop  has  an  important  influence  on  the  seismic  radiation, 
particularly  on  the  relative  excitation  of  high-  versus  low- 
frequency  signals.  This  should  be  pursued  in  the  context  of 
a  deterministic  source  model,  in  which  stress  drop  and  rupture 
velocity  are  coupled  through  a  failure  mechanism. 


VI.  DISCRIMINATION  EXPERIMENT 


Our  objective  in  the  discrimination  experiment  is  to 
analyze  short-period  seismic  waveforms  from  a  large  population 
of  events  in  order  to  identify  the  events  as  either  earth¬ 
quakes  or  underground  explosions.  During  this  reporting 
er'i0 d~we~  processed  P-wave  seismograms  fof-  52  events  recorded 
at  the  the  classified  stations. 

The  MARS  computer  program  that  is  being  used  in  the 
discrimination  experiment  is  described  in  Section  III.  During 
this  reporting  period,  two  modifications  were  made  to  the 
procedure  used  for  estimating  the  variable  frequency  magnitudes 
One  of  these  modifications  involved  the  development  of  an 
algorithm  for  estimating  the  effect  of  "local  seismic  noise" 
on  a  transient  signal  (Masso,  et  al . ,  1978).  The  second 
modification  consists  of  averaging  weighted  magnitude  estimates 
at  several  different  frequencies  over  a  low  (e.g.,  0.6  to 
0.9  Hz)  and  a  high  {e.g.,  2.5  to  3.5  Hz)  frequency  band.  The 
weights  are  based  on  the  ratio  of  signal  power  to  noise  power. 
These  weighted  mean  magnitudes  are  more  stable  estimates  than 
values  computed  at  the  individual  filter  frequencies  (Savino, 
et  al. ,  1978)  . 

In  the  following  we  will  give  a  brief  description  of  the 
changes  to  MARS  and  conclude  with  a  narrative  summary  of  the 
discrimination  results  obtained  to  date  for  the  classified  data 
set. 


6.1 


SEISMIC  NOISE  CORRECTION 


Earlier  applications  of  the  VFM  approach  to  discrimina¬ 
tion  included  a  rather  crude  noise  correction  to  the  m^ff)  data. 
The  particular  form  of  this  correction,  namely  the  subtraction 
of  a  frequency  dependent  average  noise  level,  was  observed  to 
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result  in  a  biased  (high)  estimate  of  the  probable  noise  level 
occurring  during  a  signal  time  window.  In  addition,  the  phase 
of  the  noise  relative  to  that  of  the  signal  was  not  taken  into 
account.  We  have  now  developed,  and  are  routinely  using,  a 
more  accurate  noise  correction  that  takes  advantage  of  the 
narrow-band  filtering  procedure  (Masso,  et  al.,  1978). 

In  principle,  the  noise  correction,  as  now  applied, 
consists  of  a  deterministic  component  and  a  statistical  com¬ 
ponent.  The  deterministic  component  is  based  on  a  superposed 
pulse  model  for  the  noise,  where  noise  is  defined  as  all 
energy  or  group  arrivals  not  identified  with  the  particular 
signal  being  considered.  Given  this  definition,  the  "noise" 
can  be  made  up  of  what  would  ordinarily  be  considered  signal 
(e.g.,  the  coda  of  the  first  arrival  P-wave) ,  as  well  as 
normal  background  seismic  noise.  The  particular  form  of  the 
correction  that  treats  this  "local  seismic  noise"  is  termed 
deterministic  because  the  effects  of  both  the  amplitude  and 
instantaneous  phase  of  the  "noise"  on  the  signal  can  be  cal¬ 
culated,  at  least  to  first  order. 

As  described  in  Section  III,  the  output  of  the  narrow 
band  filtering  process  consists  of  maxima  of  envelope  func¬ 
tions  as  a  function  time.  The  deterministic  noise  correction 

is  formulated  as  follows.  Let  A*(f)  be  the  measured  envelope 

g 

amplitude  associated  with  a  signal  of  interest,  and  t*(f)  be 

y 

the  energy,  or  group,  arrival  time.  In  addition,  let  (A  (f) } 

N  ^ 

be  a  set  of  noise  peaks  with  group  times  {tg(f)}  such  that, 
either:  t*  -  5t  <  t^(f)  <  t*  -  it,  or  t*  +  it  <  t^(f)  <  t*  +  5t 

g  -g-g  g  -g-g 

where : 


10  2n2  „  _  1 

Aw  ;  2Aaj 


and  suppose  that  there  are  M  such  noise  peaks.  Here  Aw  is  the 
half  power  band  width  of  the  Gaussian  filters  that  are  used  in 
MARS  (see  Section  III). 
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We  compute  the  "local  noise  corrected  signal  amplitude, 

A**(f),  from: 
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** 


(f) 


Ag  cos  fp(tg) 


exp 


•^2  2 
48  (tg  wg; 


*  it  if 

A  sin  v  (t  ) 

rr  v  ct' 


cos  ?^a) (t*) 


j]  A^m)  exp 
m-1 


•Acu  < .  *  (m)  >  2 

48  (tg  "  wg  } 


sin  (tj 


•P 


1/2 


where: 


3 


in  2 
2 


*  * 

■a  ft  i 
*P  1  g' 


*  * 

A  ( t  ) 
9  9 


(rn )  * 

■?  k  g; 


(m)  /  *  \ 

•  p  V  -g; 


=  Instantaneous  phase  of  the  signal 
at  the  envelope  peak  time  (croup 
time)  t* (f) . 

=  Amplitude  of  the  envelope  at  the 

*  .«  •  r*  "i ■  ,  — •*  '  ' 

peak  in  the  envelope  function,  occur- 
★ 

ring  at  t  ( f ) . 
g 


=  Instantaneous  phase  of  the  n  "noise 
pulse  at  the  envelope  peak  time 
(croup  time) ,  t^(f) . 


*  jj  ( t  - 
o  g 

phase  at 


'S 

(t“‘)  =  the  noise 

-  g 

the  sicr.al  grouo  time. 
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(f)  =  Amplitude  at  the  envelope  peak  for 

•»  X- 

the  m  noise  pulse  (occurring  at 
the  time  t” (f ) ) . 

Here: 


exp 


-^2 
46  1  "g 


cos 


is  the  deterministic  noise  correction  for  the 

★ 

signal  with  group  time  t^,  while 


2l 

i 

J 


si 


4m)  (O 

?  <3 


is  the  quadrature  component  corresponding  to 
We  then  compute  both  of  these  and  evaluate  the 
(deterministic)  noise  correction  I^A^I  as: 


and  save  for  later  use  in  describing  the  noise 
oopulation  in  the  m^f)  .  ane.  This  correction/ 
while  deterministic,  is  only  to  first  order. 
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<* 


** 


Associated  with  A^  (f)  will  be  an  uncertainty  due  to 

the  impossibility  of  resolving  and  correcting  for  noise  peaks 

1 


within  the  range  At  = 


2Au) 


on  either  side  of  the  signal  group 


time  (t  ),  at  a  frequency  f.  In  addition,  we  will  want  to 
9  ** 

include  in  the  uncertainty  attached  to  the  estimate  A  ,  the 

impossibility  of  making  an  exact  deterministic  correction  for 

noise  contamination  from  the  "local''  noise  pulses.  This 

uncertainty  will  be  +  AA^  with  AAj^  given  by: 

^  (T"}  AuT  (T~)  ^ 

o  o 


where 


1  =  Total  (standard)  time  window  being  pro 

cessed  (typically  above  100  sec) . 

L  =  Number  of  noise  ceaks  in  7_ .  (Use  the 
■  - 1  ■•"*--  *  o 

previously  analyzed  time  window  where 
ail  noise  peaks  have  been  identified.) 


A^(f)  =  Mean  of  the  envelope  amplitudes  at 
the  noise  peaks  in  the  window  Tq , 
i.e.  ,  A^(f)  =  £  A^Z)  (f)/L. 

1=1 


(Note  that  L  a  Aw  T  as  Au  -*■  0 ,  so  that 

o 

AA^  -*■  ^  as  Aw  -*•  0.  Further,  <5A^  -  0  as 
Aw  -*■  0  by  definition  of  the  time  interval  for 
the  5A^  correction.  Therefore  as  Aid  +  0,  and 
the  filter  Q  become  infinite,  we  get  the  usual 
noise  correction  of  a  Fourier  spectrum.)  Thus 
the  signal  spectral  amplitude  will  be  described 
by  A**(f)  +  AA^(f)  . 

89 


J 


r 


i 

? 

►v 
►* . 


k' 

k- 


»- 


The  final  step  in  this  formulation  is  the  computation 
of  the  noise  corrected  instantaneous  phase.  This  is  given  by: 

A* 

y  -  o 

y*  ~ 


** 


(t  (f))  *+tan 


-1 


where  5 


*N 


and  <5 


are  as  defined  above  while 


y  =  Ag  cos  ¥p(tg) 

A  *  *  *  * 

y  =  +  Ag  sin  fp(tg) 


where 


A  * 

y 


is  the  quadrature  signal. 


6.2  NOISE  CORRECTED  MEAN  WEIGHTED  mb(f)  ESTIMATES 

The  application  of  the  noise  correction  derived  above 
is  based  on  a  magnitude  relationship  similar  to  the  one 
originally  proposed  by  Gutenberg  and  Richter  (1956) : 

mb(f)  =  log1Q  [Af]  +  b 

where  b  is  the  distance  correction  factor.  Using  this  rela¬ 
tionship  we  compute  mt(f)  values  for  the  signal  of  interest 
from  A*  (f )  /  Ag *  ( f )  ,  Ag*(f)  +  AAn  and  Ag  (f)  -  AA^j.  The  mb(f) 
are  computed  at  frequencies  corresponding  to  the  center  fre¬ 
quencies  of  the  entire  set  of  narrow  band  filters  being  used 
(i.e.,  typically  30  filters  covering  the  band  0.5  to  5  Hz). 
For  discrimination  purposes,  however,  two  sub-bands  are 
defined:  a  low  frequency  set  { f T } ,  where  0.5  <  f  <  1.0  Hz;  a 

high  frequency  set  {f H>,  where  2.5  <_  f  £  3.5  Hz.  In  order  to 
obtain  more  stable  magnitude  estimates  than  those  based  on 
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individual  frequencies,  we  are  now  computing  mean  weighted 
m^ff)  values  over  the  frequency  sets  {f^}  and  {fH>  and  using 
these  in  the  discrimination  experiment.  The  magnitude  rela¬ 
tionship  is  given  by: 


A* 

V£)  *  losio  I  § w*  (ATfk> '  ^  wk 


Vi 

E  w„  ♦ 


where 


„k(£)  =  [A*/A„]2 


are  weight  factors,  measuring  the  signal  information  content 
or  "quality"  of  the  signal  information  at  the  frequency  f^. 
The  summation  is  performed  over  the  low,  and  high 

nVj_> ( ffl)  ,  frequency  discrimination  sub-bands. 

Associated  with  the  magnitudes  m.(fL)  and  ^(fjj)  are 
the  frequencies  and  fH  defined  by: 


w  N 

iT  =  E  w,  f^k)  7Ewi 

L  k=l  K  k=l  * 


IN 

-  2  wk  4k)  7  IE  ^ 


Finally,  we  also  compute  the  uncertainty  in  the 
weighted  mean  magnitudes  at  and  fH  due  to  the  uncertainty 
in  the  noise  correction.  These  are  given  by: 


J 


Note  that  Am^  ?£  -  Am^  in  general,  and  that  none  of  the  An^ 
values  obtained  for  the  two  sets  of  discrimination  fre¬ 
quencies  will  be  equal  in  general.  We  get  from  this  then 
four  distinct  An^  values. 

The  results  of  applying  the  newly  formulated  noise 
correction  to  event  m^ff)  data  are  summarized  in  Figures  25a 
and  25b.  Figure  25a  shows  the  behavior  of  typical  event 
(uncorrected  for  noise)  and  noise  populations  in  the  m^ff) 
plane.  This  figure  is  a  generalization  of  earlier  results 
previously  reported  on  (Savino,  et  al. ,  1975;  Rodi,  et  al. , 
1978) .  Figure  25b  demonstrates  the  manner  in  which  the 
earthquake  and  explosion  populations  separate  when  the 
deterministic  and  statistical  noise  corrections,  together 
with  the  uncertainty  inherent  in  both  these  corrections,  are 
applied.  The  enhanced  separation  of  populations  is  especially 
significant  in  the  low  m^f)  range  where  noise  plays  an 
important  role.  The  primed  mb(f)  values  in  Figure  25b  refer 
to  weighted  mean  magnitudes  determined  from  the  low  and  high 
frequency  bands.  This  is  the  procedure  that  we  are  routinely 
using  in  the  discrimination  experiment  for  the  computation  of 
the  variable  frequency  magnitudes. 
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LOW  FREQUENCY  BODY  WAVE  MAGNITUDE  [m/(f_)l 


Figure  25a.  Typical  event  distributions  in  the  itivj  (f)  plane  for 

event  data  that  is  not  corrected  for  noise  contamina¬ 
tion.  Noise  pulses,  when  viewed  in  this  space  appear 
roughly  as  shewn  and  affect  explosion  event  ( f  ^ ) 

values  most  strongly,  causing  population  overlap^at 
low  magnitudes.  The  population  boundaries  for  noise 
and  events  are  somewhat  source  and  receiver  dependent 
due  to  earth  structure  variations. 


» 
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FREQUENCY  BODY  WAVE  MAGNITUDE  [m,*(fT)] 


c 


ICN  B 


HIGH  FREQUENCY  BODY  WAVE  MAGNITUDE  [m£(fH) 


Figure  25b. 


Typical  event  distributions  in  the  mv(f)  plane  for 
noise  corrected  event  data.  A  discrimination  line 
can  be  defined  on  the  basis  of  the  definition  of 
these  populations  using  known  events,  or  on  the 
basis  of  theoretical  predictions. 


6.3 


DISCRIMINATION  RESULTS 


Event  seismograms  for  52  events  recorded  at  one  or  more 
of  the  eight  classified  stations  were  received  at  S3  and  pro¬ 
cessed  during  this  reporting  period.  Preliminary  results  indi¬ 
cate  that  the  variable  frequency  magnitude  (VFM)  approach  can 
discriminate  events  down  to  small  magnitude  levels  (low  P-wave 
signal- to-noise  ratios)  at  those  stations  characterized  by 
relatively  low  background  noise  levels.  There  is  also  a 
suggestion  that  stations  located  over  probable  high-Q  upper 
mantle  regions  (e.g.,  shields)  provide  better  separation  of 
earthquakes  and  explosions  than  do  stations  located  in  tec¬ 
tonic  regions  of  high  heat  flow  and  large  positive  (slow) 
travel-time  residuals  (by  inference,  low-Q  upper  mantle 
regions) .  These  tentative  results  will  be  examined  in  much 
greater  detail  once  the  entire  data  base  (both  classified  and 
unclassified  stations)  has  been  analyzed. 

The  numbers  of  events  presently  at  S3  are  92  for  the 
classified  stations  and  117  for  the  unclassified  stations. 

Our  future  plans  are  to  first  complete  analysis  of  the 
classified  data  and  then  proceed  to  the  unclassified  set. 

The  data  processing  (i.e.,  MARS  runs  on  all  the  seismograms) 
should  be  completed  by  the  end  of  February.  Assuming  a 
modest  number  of  additional  events,  final  reports  on  the 
Priority  1  and  2  station  sets  will  be  completed  by  mid- 
April.  _ _  _ 
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